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INTRODUCTION  AND  STATEMENT  OF  PURPOSE 

New  methods  for  design  and  "i&nagement  of  water- resource 
systsns  are  being  evolved  as  part  of  a  general  social  tendency 
toward  expressing  social  problems  in  the  formal  models  that  have 
hitherto  been  restricted  to  scientific  and  engineering  problems. 

Two  general  types  of  models  have  been  fruitful  in  the  field 
of  water-resource  development!  the  simulation  model  and  the 
analytic  model.   Simulation  is  used  for  highly  complex  systems 
involving  a  large  number  of  design  decisions.   For  less 
complicated  systems,  various  optimization  techniques  can  be 
used  to  obtain  the  optimal  conditions.   Simulations  are  awkward 
when  a  wide  range  of  design  decisions  has  to  be  evaluated: 
analytic  models  can  not  be  applied  to  practical  problems  without 
drastically  simplifying'  them.   But  the  two  nethods  can  be  used 
in  tandem,  with  analytic  models  delimiting  the  range  within 
which  simulation  is  required. 

Two  basic  problems  in  the  analysis  of  water-resource 
systems  arei  the  problem  of  obtaining  a  realistic  mathematical 

1  and  the  problem  of  obtaining  the  best  management  policy 
from  an  analysis  of  the  resulting  aathematioal  model. 

Of  the  different  analytical  approaches,  dynamic  programming 
has  been  quite  useful  In  optimization  of  water-resource  systems. 
Linear  nrograamlng  is  not  generally  useful,  bees-use  individual 
reservoir  utility  or  cost  functions  are  in  general  definitely 


considered . 

Significant  progress  has  been  r 
nonlinear  programming  problem  durinj 
different  approaches  are  generally  i 


de  in  treating  the  gener.' 
the  last  few  years.  The 
combination  of  gradient  < 


search  techniques  with  a 

specific  purpose  of  this 
optimization  technique  - 


lod  for  treating  const: 
sarch  is  to  apply  a  rei 
:n '  s  a rad  i  ent  pro i  ecti < 


ints.   The 
ntly  developed 
method  [2^]  - 


to  the  optimal  management  and  design  of  water-resource  systems. 

The  gradient  projection  method  is  described  in  Chapter  2. 
A  rather  nontheoretical  approach  is  attempted  in  describing 
various  aspects  of  the  method.   The  emphasis  is  on  the  applicatlio: 


of  the  method  and  nc 
Two  test  systet 


ir.h'.n  the  framework  of 
The  basic  principle 


the  complex  mathematical  details, 
me  connected  with  water-quality  and 
!r-quantlty  -  are  proposed  and  solved 
m  analytical  approach, 
involved  in  the  construction  of  a 


mathematical  model  of  dissolved  oxyge 


simplified  rlv. 


are  presented  in  Chapter  3.  . 
how  the  mathematical  model  cai 
quality  control  information. 

Chapter  ^  deals  with  wat< 
fled  river  basin  cor  "iguratioi 
mathematical  progr 

'    at   this  :  ..  atem  la 


i  attempt  has  been  made  to  indicate 
be  used  to  generate  useful  water 

quantity  aspects.   A  very  sirapll- 
ls  selected  to  illustrate  the 
ch.   An  optimal  plan  for 


THE  GRADIENT  PROJECTION  METHOD 

2.1    INTRODUCTION 

In  the  fields  of  nonlinear  optimization  there  are  many 
methods  [29]  which  may  be  used  successfully  on  unconstrained 
problems,   However,  most  practical  problems  involve  constraints, 
both  general  and  specific,  which  must  be  satisfied  in  order  to 
obtain  meaningful  results. 

Rosen  [24]  developed  the  gradient  projection  (OP)  method  for 
solving  the  subclass  of  nonlinear  programming  problems 
characterized  by  linear  constraints.   He  extended  this  work  to 
handle  the  more  difficult  case  of  nonlinear  constraints  [15].   The 
discussion  here  is  confined  to  the  simpler  case  of  linear 
constraints  only. 

Briefly,  GP  is  a  method  for  a  system  of  m  variables  (x's) 
subject  to  linear  constraints  -which  may  consist  of  inequalities, 
equalities  or  both.   The  objective  is  to  find  the  maximum  value 
of  the  objective  function  while  satisfying  the  constraints. 
In  geometric  terms,  GP  starts  with  a.  feasible  point  (one  that 
fritisfies  all  the  constraints)  or  finds  one,  if  not  given,  and 
then  a  ntep-'ise  procedure  gives  a  new  feasible  point  at  each 
•  top  with  an  Increased  value  of  the  objective  function.   This  is 
accomplish:-:  "...  ■  -'  -.-.-  ..-.ps  in  the  direction  of  the  gradient  of 
the  objective  function  or  its  projection  on  the  Interact '  r.i  of 
Selected  constraints.  Th  i  maximi   ;  -   -   i       -..-   ;ep  that 
would  increase  ths  r;-:lus  of  the  obj&ctive  "unction  would  violate 


the  constraints,  or  when  the  gradient  roes  to  zero. 

2.2   .-OIt"uIA.riCN  FO:?  LTKSA3  CC'TSTP.AIIJTS  [24-] 

A   problem  with  m  variables,  x ,  i  =  1,  2,  ...,  a,  is 
considered.   Geometrically,  any  specified  set  of  values  for  the 
x.     represents  a  point  in  a  Eu.elld.ean  ™~  dim  ens'  onal  erace,  E  . 
It  5  s  assume-]  that  the  variables  are  constrained  by  a  set  of  k 
linear  inequalities  or  equalities.   These  constraints  form  a 
convex  region  ?.  in  E„  which  is  assuued  to  be  bounded.   A  feasible 
point  lies  in  S.   The  constraints  are  of  the  form 

jx  vrbiV         l"1' 2 k  (1) 

£      (n=  J2  -.   1  i  «  1,    2,    ....   It  (2) 


3b.ce   a  is   a  bounded   region  there  mat 
so  that  k  >  m*l. 

Cc  >OMding   to   each   of  the  k   c 

defined 

tti  =   Ki'\2 \*\ 


now   be  writti 


!   by  I,;-)    :01e< 


hyperplasia  whloh  Will  be  denoted  by  H  ,  m 

H  t  ^1(x)  =0         1  „  1,  2,  ....  k  (5) 

The  unit  vector  n  is  orthogonal  to  H,.   A  set  of  hyperplanes  are 
linearly  independent  if  the  corresponding  vectors  n  are  linearly 
independent.   Clearly  there  can  be  at  most  a  linearly  independent 
hyperf planes  in  an  m-dimensional  space. 
Define  the  m  x  k  matrix 

and  the  k-dlmensional  vector 

\-     [\>    \ \\  (7) 

Then  (.1)  or  {h)    oan  be  written  conveniently  as 

N*  x  -  \  >   °  (3) 

A  set  of  q  linearly  in~i      t  unit  v  ctors  n  . 
1 
1  =  li  2 q,  with  1  <  q  <  a,  defines  a  set  of  q    linearly 

independent  hyperplanes  1^  as  given  by  (4)  and  (5)  for  any 

specified  iralBes  of  Wie  b  .  Let 

VCVB2 V  (9) 

It  cct  be  shown  that  the  a   x  q  symmetric  matrix  tJ3   H     is 

■'■-■-■  ,"10  j    r-r.d    thei-,:-foi  e    it:;    ':iv;r,;=    fi;T   \i    } -1    ,, -fist's 

:.'.■  intersection  of  the  q   hyperplanes  H     and    let 
~i  be  the  q-dlmenslonal  subspacs  of  E      ihl    h   is   spanned   by 


;  "space.   The  suhspaces  Q  and.  Q  sre  orthogonal. 
Define  the  m  x  a   symmetric  aatrix 


i  project!. on  'natr: 


An  a   ::  r.  L2:<trix   is   now  defined  by 

Pq  =   I   -  Pq  (U) 

The  matrix  P   is  a  projection  nstrix  which  takes  any  vector  in 

3  into  the  intersection  Q. 

In  course  of  an  optimization  calculation,  it  is  necessary 

to  obtain  the  projection  of  the  gradient  vector  on  various 

intersections  Q.  In  other  words,  the  matrix  (A;-1  la  required 

at  each  step.   Rosen  [24-]  derived  two  recursion  relations 

which  permit  a  hyperplane  to  be  dropped  from  or  added  to  (N  N)"1 

with  considerably  lews  computation.  A  simple  and  verj  useful 

recursion  relation  for  ?   is  also  given. 
1 

2.3   CO!  .  ;  [iV] 

•  . active  function 

■•:■;  -  '(*!•  ->:z V  (12) 

is    defined   in  R,   and   is  \  \  ■  -:-A 

ipect  to    .he  x.  .      Msi       [x) 
■  find  a  La 


B  at  which  F(x)  has  a  global  maximum.   If  the  point  x    is  on 
the  boundary  of  5,  it  is  a  constrained  maximum. 

The  gradient  of  F(x)  Kill  be  denoted  by  g<x)  and  is  defined 
by 

s(x)  .  lL,  JS. AH  (13) 

It  is  sell  kr.rx.i  that  for  an  unconstrained  concave  function 
the  necessary  and  sufficient  condition  that  x  be  the  maximum  point 
is  g(xo)  a  0.   If  such  a  point  exists  in  the  interior  of  P.,  it 
Will  be  called  the  interior  global  maximum  of  P(x).   If  g(x) 
does  not  vanish  in  the  interior  of  R,  the  global  maximum  lies  on 
the  boundary  and  will  be  called  a  constrained  global  raximura. 
The  basic  theorem  concerning  a  constrained  global  maximum  is  as 
follows 
Theorem i 

Let  x  be  a  boundary  point  of  B  which  lies  on  exactly  q, 
1  <   q  <  m,  hyperplanes,  which  are  assumed  to  be  linearly 
independent.   Let  the  intersection  of  these  hyperplanes  be  the 
manifold  Q.  Then  the  poiat  x  is  a  constrained  global  maximum 
of  F(x)  If,  and  only  if, 

P  g(xj  -,   0  (14) 


«<x J  <   0 


I  n  a    th,  on  a    to     st    jllsl      he 
snce  to  a   global  maxlraui    of  F(x).      Poi 


Ta-„[-.c:Eatic<3l  r.j-fli.tr.er.t  of  these  ^spool's  of  constrained  -aximusE 
K.-y-  c-nv-r  •.  .ca,  the  interested  reader  csn  consult  the  original 
reference  [»]. 

CI  [M  [24] 
Certain  quantities  used  in  the  algorithm  are  defined  here 
Without  going  into  complex  mathematical  details.   Let 

%  -  *<V  (16) 

to  a  constraint 

actor  in  direction  of  step     (18) 

Y     =  selected  stec  size 

i:o\-  consider  a  feasible  point  x  which  lies  on  the  a 
0  ' 

\.  i  =  1 q.   Since  it  is  a  feasible 

point,  X.1(x)  >  0,  1  =  q+l,  ....  k.  Then  for  each  of  the  remaining 
k-q  hyperplwies,  iy  1  „  q.+l ':.  there  iay  exist  a  value 

V  =  'i,  such  t:a:  >.  (.■:)  a  Q.     ?    ls  fche  distance  from  x  to  the 

0 
fallal  to  2.   In  particular, 

r,  ■  J~       1  --  <H-1,  ....  "  (19) 

i  •  .  c  of  r, 


urn 

sal   distani 

Lri 

'    r2 

^r 

S0           ' 

£cli 

Tr 


VK^ol 


"1  =   "0  T  '  m' 


ra   ■    =la|Yi    >    0]  i   =    q~l fc  (20) 

The  distance  Ta  represents  the  largest  step  that  can  be  taken 
from  y_0  in  the  direction  2  without  leaving  H  (see  Fig.  1).   if 

*  =  -o  +r*  (2D 

then  it  follows  that  for  any  Y,   0  <  f  <  f  i  X  is  in  E. 

The  algorithm  stated  below  is  for  a  current,  arbltary  point 
x  in  H  which  lies  on  the  manifold  Q  forced  by  the  intersection 

of  q  linearly  independent  hyperplanes. 

Step  1. ' 

Compute  s^  =  a-(*v)  and  P^.   If  p  g?  „  o  and  r  <  0,  where 
r  is  given  by  (17),  then  x^  is  a  global  mftxim.ua. 

Step  2. 

,oW>ieither||Vv|>_jo.lr     dqi--l)orj|PqM< 

1  ~2  -5        -i 

2  rq  'qq   >  Khere  rq  dqq  2  >   P4  &Xi    2 ,  1  =--  1 ,  2 q-1,  and 

where  a   is  the,  1th  a  1  agonal  element  of  (NT?J  f\   rf  the 
11  q  q 

■  lolds,  compute  z  according  to  (18).  If  the  latter  holds, 

'   '" '  ■■  :   °  '' '  ~  Vi  ■  '    '   eo*  '■  '  "  :'  F,-i- 

Compute 

z    .      ?: 

'  '  K^  1 

Strip    3. 

V  .     'rora    (19)    Wd    (20) 

Sto  1    ■'. 


Compute   g*        =   g(x'       ).      If   zTg'    ,    >   C,    take  x  =   x' 

v+1  v+1  v+1  -  v+l  v+1 

and  add   to  Q  the  H.    corresponding   to  the  minimum'1'     in   (20). 


Vl*  ?K+1+    (1  "f)xv     •  (23) 

The  intersection  Q  remains  unchanged. 

This  is  the  basic  algorithm  for  a  single  step  from  a  point 

x   (if  it  is  not  maximum)  to  a  new  point  x    with  an  increased 

v  v+1 

value  of  Fix).      At  each  step,  the  intersection  Q  may  remain 
unchanged,  a  hyperplane  nay  be  added  or  dropped,  or  one   may  be 
dropped  and  another  added.   The  last  step  of  algorithm  is  the 
orthogonal  gradient  interpolation. 

2.5   PROGRAM  FLOW 

A  genera]  ■.  .     .  .   ■  ,  known  as    ■■  c  [  >l]  was  developed 
for  the  purpose  of  oclving  a  ser;  •■  ■•  •  -     .   ,--rig  design  and 
,;';c;''""'"  i:r.-!c3cr.s.   A  complete   loi  chart  -.'  the  nrojram  is 

pendlx  B.   This  .'.  -  '  -elation  methods, 

matrix  computations,  re-lnversj  i  lo  Lo  jto. 

A  brief  description  of  th«  md  subroutines 

follows  the  explanation  of  I.  ol/;r  .-;:■:-.-•.=  .?nd  li.iits  Meces^.ry  for 


the  program. 

2.5.1   Tolerances  and  Limits 

The  gradient  tolerance,  g^  is  used  to  determine  when  the 
norm  of  the  gradient  is  zero,  and  the  problem  has  reached  a 
maximum.   The  value  of  £1  is  harder  to  determine  for  a  nonlinear 
function.   In  general,  the  smaller  the  value  of  t    ,  the  better 
the  "maximum"  will  be.   The  price  paid  for  this  better  answer 
will  be  in  more  machine  tine.   A  value  of  I      about  10~3  ||r||  seems 
to  be  reasonable. 

The  constraint  tolerance,  £2 ,  determines  when  a  point  is  on 
a  constraint,  and  is  therefore  the  acceptable  error  in  satisfying 
the  constraints.   A  value  of  £  about  10"3  V ,  where  b'  is  the 
largest  right  hand  side  of  constraints,  b  ,  divided  by  the 
corresponding  scaling  divisor  or  10-3  times  the  largest  value  of 
x  seems  to  be  reasonable. 

The  linear  dependence  tolerance,  £,,  determines  when  a 
constraint  is  linearly  dependent  and  therefore  can  not  be  added 
to  the  basis.   The  general  value  used  for  £  is  0.005. 

^max  ls  the  maximual  step  length  which  is  used  as  the  initial 
step  length  and  must  be  provided  when  the  region  is  unbounded. 
A  possible  choice  for  7^  Is  IJn,    where  L   is  the  largest  value 

that  srvj   Individual  s  oan  assume.  A  value  off    -  10  1  - 

•max 
found  reasonable  for  quadratic  functions.   The  program  computes 
a  minimum  step  length,  f         „  lo"'1'  f       ,  which  prevents  some 
Ing  or  extrapolating  when  the  distance  between  the 

points  is  less  than  T 

1  win 

7"  :  five  111  its  which  must  be  set  are  now  su-unarized.   The 
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settings  for  these  limits  will  generally  depend  on  the  size  and 
type  of  problem  being  solved. 

MAXU,  maximum  number  of  steps.   The  number  of  steps  required 
to  reach  the  maximum  is  difficult  to  pre-determine  since  it 
depends  on  such  factors  as  size,  type  of  function,  and  number  of 
constraints  in  the  basis.   It  keeps  the  run  to  a  reasonable 
length  without  having  to  rely  on  a  time  limit. 

MXBNi  maximum  number  of  re-inversions.  The  number  of  re- 
inversions  is  limited  to  prevent  the  problem  from  re-inverting 
too  often  thereby  consuming  extra  time.  The  program  will  not  re- 
invert  twice  in  a  row  to  the  same  basis,  but  it  may  re-invert 
after  only  one  basis  change  or  possibly  repeat  a  series  of  re- 
inversions. 

^raax  and  Yrnax  ■1-lmit  the  n™ber  of  gradient  interpolations 
which  are  computed  in  finding  an  x  for  which  F  has  increased. 

For  a  quadr.-itic  function,  8     should  be  one  and  Y    should  be 

IEcix  max 

zero"  Tnax  Umlts  the  number  of  points  to  be  saved  and  used  in 
computing  extrapolation.  If  %3y.  Is  zero  in  input,  the  program 
will  use  t;;-i  theoretical  limit,  m-q,  fori?    . 

2.5-2   Main  Program 

"he  fl:.-.-:c  ;■■  action  reals  input  and  s.-?ts  the  initial  conditions 
for  a  pr.Volem.   Tile  constraints  are  read  end  normalized  to  unit 
Vectors.   If  equalities  are  todicated,  ;,he«  form  the  initial 
basis',  and  the  corresponding  inverse  is  soaputed. 

In  the  next  section,  tn,v.  urogram  t-^'.s  for  fe^siMlltj 

./,  finds  a  feasible  x  or  determines  that  there  la 


Ik 

no  feasible  x. 

When  a  feasible  x  is  found,  the  subprogram  is  entered  to 
compute  F  and  g,   The  program  is  now  ready  to  enter  the  step 
procedure.   A  step  includes  the  projection  of  the  gradient, 
testing  for  the  maximum,  changes  in  the  basis  and  computation  of 
the  step  length. 

At  the  beginning  of  each  step,  the  program  has  a  feasible 
x  for  which  P  and  g  have  been  computed.   The  non-basis  constraints 
are  classified  into  V(  (xj  <  £  g)  and  W(  |\J  >  g  ) .   The  norm  of 
the  gradient  is  computed  and  tested  for  zero.   If  flgjl  <  g  , 
this  is  the  maximum. 

If  ||'|  >  £1   and  q  =  0,  the  computed  gradient  is  stored  as 
the  projected  gradient,  and  program  skips  to  comptite  z.   If 
q  #  0,  the  ^racUent  is  projected  and  the  norm  of  the  projected 
,'jradient  is  tested  for  zero.   If  ||ps||  <  £    the  program  tests 
for  a  constraint  to  drop.   If  there  is  no  constraint  to  drop, 
the  maximum  test  is  satisfied  and  the  current  point  is  the 
maximum.   If  there  Is  at  least  one,  the  best  one  is  cropped  from 
the  basis,  and  gP   HI  is  calculated  and  tested. 

If  ||F-'|  >€   and  q  <  ra,  the  unit  vector  z  is  conputed.   If 
V  j.   0,  7.   n  1:5  computed  and  tested  for  all  i  in  V.   If  the 

IH  z  n,  is  negative,  the  corresponding  constraint  is  added 
to  the  basis  if  it  is  linearly  independent,  and  |Pg|  is  calculated 
and  tested. 

When  there  are  no  more  constraints  in  V  for  which  zTn.  is 
—  v  -  0,  the  program  tests  for  basis  ohanges. 
If  the  rtep  Mas  interior,  the  program  tests  for  a  constraint  drop. 
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When  the  program  finds  no  more  changes  to  be  made  In  the 
basis,  it  is  ready  to  compute  the  next  step  length.   If  the 
previous  step  was  not  interior,  the  initial  step  length  Is  set 
at  ^max"   If  the  ?revi0lis  steP  was  interior,  the  initial  step 
length  depends  on  the  extrapolation  method  selected.   Two  extra- 
polation methods  available  are  x,  G  method  and  3-point  method. 

When  initial  step  length  is  determined  and  if  W  =  0  this 
is  checked  against  constraints  in  W.   For  all  1  in  W  for  which 
z^n^   is  negrative,  Xt   =  -^A^  is  computed.   If  the  minimum  1^ 
is  less  than  the  initial  T,  the  step  is  limited  by  the  constraint 
H1  and  T,  replaces  Y. 

A  new  x,  x^  =  ae^  +  fz,  is  computed.  F  and  g  are  computed 
and  zTn,.  is  computed.  If  T  is  less  than  Yml  ,  the  interpolations 
axe  skipped.  Two  kinds  of  interpolations  are  available  to  find 
the  maximum  of  F  in  the  direction  determined  by  z.  The  first  is 
computed  if  the  gradient  has  reversed  direction  between  x  and 
x^.  This  prevents  the  problem  from  overshooting  the  maximum. 
The  second  is  computed  if  F  has  not  increased. 

After  completing  the  necessary  gradient  interpolations,  x  is 
checked  for  feasibility.   Theoretically,  x  should  be  feasible, 
but  because  of  inaccuracies  due  to  round-off  in  the  inverse 
matrix  it  Is  possible  to  violate  the  constraints.   Tn  that  case 
F  and  g  are  recomputed  and  ?  is  rcchecked. 

One  pass  through  the  step  procedure  is  now  complete.  The 
step  counter  is  tested  against  maximum  step  limit.  If  the 
limit  is  not  reached,  the  program  returns  to  the  beginning  of 
step  procedure. 


2.5-3   Subroutines 

Subrpitone  REINV  essentially  computes  the  inverse  matrix 

T   -1 
(N  N  )   whereas  subrouLir.-s  ::..:cc;:  and  CC;.:AT  do  the  macrix 

q  q 
computations  required  by  the  gradient  projection  algorithm. 

Subroutine  AMOA  calculates  the  \*s  while  subroutine  CLASS 

classifies  the  constraints  into  different  catagories.   In  the 

program. 

u  •-  linearly  dependent  constraint 

v  =  constraints  not  in  the  basis  with  X   =  0 

w  =  constraints  not  in  the  basis  with  X  >   0 

Subroutine  FUNGI  is  added  to  the  program  to  compute  the 
value  of  the  objective  function,  F(x),  and  its  gradient,  g(x), 
for  any  given  x  within  the  region. 

2.6   DISCUSSION 

The  technique  and  program  into  which  it  is  incorporated  are 
designed  to  handle  general  nonlinear  problems  with  linear 
constraints.  The  method  is  intended  to  apply  to  a  variety  of 
problems,  wen  though  in  many  cases  more  computation  time  may 
be  required. 

If  the  function  is  concave,  a  global  maximum  is  guaranteed. 
In  order  cases,  the  solution  may  be  only  a  local  maximum,  and 
several  widely  separated  starting  points  should  be  tried.   If 
different  results  are  obtained,  the  best  that  can  be  done  is  to 
take  the  maximum  value  of  this  set. 

As  true  with  different  variants  of  the  gradient  methods,  in 


17 

general,  an  Infinite  number  of  iterations  may  be  required  before 
the  conditions  regarding  a  constrained  global  maxinum  are 
reached.   Also,  the  steps  may  get  too  small  near  the  stationary 
point  resulting  in  a  very  slow  convergence  rate  [ll]. 

Goldfarb  [9]  has  developed  a  conjugate  gradient  method  for 
nonlinear  problems  with  linear  constraints.   According  to  him, 
the  method  has  a  faster  convergence  rate  and  can  navigate  through 
the  troublesome  regions  like  "steep  ridges"  and  "narrow  curving 
valleys".   At  this  stage,  only  a  limited  amount  of  computational 
experience  is  available  for  conjugate  gradient  method. 

In  conclusion,  with  this  program,  as  with  most  nonlinear 
programs,  the  results  obtained  on  real-life  problems  are  very 
dependent  on  the  design  of  the  problem  to  be  solved,  as  well  as 
the  effectiveness  of  minor  adjustments  of  an  algorithm  to  obtain 
best  results  for  specific  unusual  problems. 


A  MANAGEMENT  MODEL  FOR  WATER  QUALITY  CONTROL 

3.1   INTRODUCTION 

The  problem  of  river  basin  planning  for  water  quality  is 
currently  receiving  wide  attention.   The  stimulus  for  this 
attention  is  the  recognition  that  water  quality  problems  are  not 
necessarily  the  result  of  one  recalcitrant  polluter  but  of  many. 
Because  of  this,  the  Federal  government's  deep  commitment  to 
restoring  stream  quality  is  increasing  not  only  in  form  of 
financial  assistance,  but  also  in  responsibility. 

One  municipality,  one  industry  or  even  one  state  can  not 
always  control  stream  pollution,   For  this  reason  governmental 
agencies  having  jurisdiction  over  entire  ri^er  basins  are 
establishing  quality  standards  for  each  section  of  the  stream. 
These  stream  standards  are  intended  to  maintain  stream  quality  by 
limiting  the  amount  of  waste  that  can  be  discharged  into  the 
stream. 

The  problem  of  determining  standards  for  the  stream  becomes 
more  complex  when  there  are  two  or  more  sources  of  pollution. 
In  these  cases,  the  amount  of  waste  released  from  one  point  may 
mix  with  the  waste  released  at  another  point  to  contribute  to 
the  pollution  downstream  from  both  points.   The  quality  standard 
at  any  point  in  the  stream  can  be  met  by  many  combinations  of 
quantities  released  at  various  locations  upstream.   The  problem 
is  to  find  the  combination  that  results  in  a  minimum  total  cost. 
Systems  analysis  and  its  applications  have  been  increasingly 
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applied  to  this  problem  of  water  quality  management.   A  series 
linear  programming  models  were  structured  to  determine  alternative 
ways  of  meeting  quality  standards.   Deininger  [>]  structured  the 
problem  as  a  linear  programming  model  utilizing  various 
approximations  of  the  differential  equations  used  to  describe 
the  dissolved  oxygen  profile  of  streams.   Kerri  [17]  proposed  a 
dynamic  model  for  achieving  and  maintaining  water  quality  control. 
He  applied  the  concept  of  a  critical  reach  to  a  simplified 
version  of  Willamette  River  in  Oregon.   Thomann  [28]  and  Sobel 
[26]  also  developed  linear  programming  models  for  related  but 
different  conditions  of  the  Delaware  estury.   Liebman  and  Lynn 
[18]  presented,  a  dynamic  programming  model  to  study  this  problem. 
The  model  was  solved  for  a  simplified  example  based  on  data  from 
the  Willamette  River.   Revelle,  Louclcs  and  Lynn  [19,  23]  developed 
several  linear  programming  models.   The  difference  in  each  model 
reflected  both  the  assumption  regarding  the  river  basin  and  the 
manner  in  which  the  standard  is  specified. 

In  general,  the  treatment  plant  costs  are  not  linear.   An 
attempt  has  been  made  to  formulate  this  problem  as  a  multistage 
decision  process  with  a  nonlinear  cost  function.   The  model  so 
developed  can  be  readily  applied  to  a  variety  of  river  basins 
with  minimal  alteration  of  the  basic  model.   The  treatment 
closely  follows  one  discussed  in  [19]. 

3-2    BACKGROUND 

A  large  portion  of  the  municipal  and  industrial  waste 
released  into  streams  is  organic  material.   These  organics 


become  a  source  of  nutrients  for  many  organisms  found  in  streams, 
Dissolved  oxygen  (TO)  contained  is  withdrawn  by  these  organisms 
in  the  process  of  utilizing  these  wastes.   The  larger  the 
quantity  of  these  bio-degradable  wastes  the  larger  is  the 
population  of  those  organisms  and  the  greater  is  the  demand  for 
oxygen. 

Fish  and  other  aquatic  animals  and  plants  require  certain 
minimum  concentrations  of  DO  if  they  are  to  survive  in  the 
stream.   If  the  DO  concentration  is  completely  depleted,  the 
stream  becomes  anaerobic  and  may  be  more  reminiscent  of  a  sewer 
than  a  stream.   Insufficient  removal  of  organics  in  the  waste- 
water prior  to  its  release  into  the  stream  can  bring  about  these 
conditions.   For  these  reasons  stream  quality  standards  usually 
specify  minimum  allowable  DO  concentrations  in  each  section  of 
the  stream.   The  CO  parameter  Is  the  one  most  commonly  used  to 
measure  and  limit  the  amount  of  pollution  resulting  from  these 
organic  wastes. 

The  capacity  of  the  stream  to  assimilate  bio-degradable 
wastes  is  determined  by  such  factors  as  stream  flow,  stream 
temperature,  the  waste  concentration  as  measured  by  its  bio- 
chemical oxygen  demand  (BOD),  the  DO  concentration,  and  the 
physical  and  biological  properties  of  the  stream  that  affect 
settling  rates,  reaeration,  BOD  addition  due  to  runoff,  scour  etc 
Two  differential  equations  describing  the  processes  that 
determine  the  concentration  of  DO  in  the  stream  have  been 
developed  by  Camp  [2]  and  Dobbins  [5]. 

The  first  equation  assumes  that  the  rate  of  change  in  the 


BOD  concentration  with  tla»,  (dB)/(dt),  is  proportional  to  the 
concentration  of  HOD  present,  3,  and  to  the  rate  of  BOD  addition, 
R,  due  to  runoff  and  scour. 

(d3)/(dt)  =  -  (k  +  k  )  B  +  R  {1) 

The  terms  J^  and  k^  represent  rate  constants  for  deoxygenatlon 
and  sedimentation,  respectively. 

The  second  equation  assumes  that  the  rate  of  change  in  DO 
deficit,  (dD)/(dt),  is  proportional  to  the  concentration  of  30D 
present,  B|  the  existing  oxygen  deficit,  D;  and  the  rate  of  oxygen 
production  or  reduction,  A,  due  to  plant  photosynthesis  and 
respiration. 

(dD)/(dt)  i   kjB  -  k2D  -  A  (2) 

The  constant,  k.,,  reflects  the  rate  at  which  BO  is  returned  to 
the  stream  through  reaeration. 

The  equations,  at  best,  grossly  describe  the  effects  of  the 
introduction  of  unstable  oxygen-demanding  substances  upon  the 
oxygen  resources  of  the  stream.   They  do  not  adequately  describe 
the  complex  biological,  physical  and  chemical  phenomena  of 
streams.   But  the  formulations  are  currently  used  by  state  and 
federal  officials  to  prescribe  levels  of  wastewater  treatment. 

By  Integrating  equation  (1),  the  30D  concentration,  B  , 
at  any  point  (corresponding  to  a  time,  t)  downstream  from  an 
initial  BOD  concentration,  3o>  can  be  determined. 


Using  equation  (J),  eq' 
the  oxygen  deficit,  D 
initial  oxygen  deficit 


At  any  time,  t,  the  DO  saturation  concei 
deficit,  D  ,  yields  the  DO  concentratloi 


"oxygen  sag"  curve  as  shown  in 
,  and  critical  fclke,  t  ,  occur 
cv;est  value,  C^.      The  critical 
;  between  saturation  concentration, 
1  at  the  critical  tine.   I-n 
.on  exceeds  the  reaeration  rate. 
In  Region  II  the  reverse  is  true. 

It  is  recognised  that  methods  for  measuring  some  of  pa.ramete 
namely  A,  k,  and  R,  have  not  been  perfected  and  in  most  cases 
these  are  unavailable.   If  a.,  k  and  R  are  assumed  to  be  zero, 
the  more  general  o:-y;-,m   sag  sqatibn  (4)  becomes  the  simpler 
Streeter-Phelps  [27]  sag  equation. 
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Fie.    2.      Oxygen   S<xZ  Curve 


3.3    MATHEMATICAL  70HKULATI0N 

The  purpose  of  this  sectior 
Linear  oxygen  sag  equation  {k)  c 
linear  constraints.   These  lines 


;o  illustrate  how  thi 
:  recast  into  a  serii 
straints  c 


be 


icorporated  into  £ 

ist  solutions  of  \ 


lathematical  model  for  determining  minimum 
■ious  water  quality  control  poticies  in  a 


A  stream  in  which  there  are  N  waste  dischanges  is  considered 

in  the  model  as  being  divided  into  N  reaches,  the  reach 

being  defined  as  the  stretch  of  stream  between  i   and 

(i+l)st  discharge. 

Tributaries,  if  any  are  assumed  to  enter  the  top  cf  a  reach. 

The  amount  of  flow  from  each  discharge,  the  DO  concentrations, 

and  the  raw  30D  loadings  are  assumed  to  be  known. 

The  parameters  of  the  equation  relating  DO  to  waste  loading 


are  constant  In  each  reach  and  known. 

5.  The  stream  flow  in  each  reach  is  considered  deterministic 
and  known. 

6.  The  standards  for  minimum  DO  concentration  In  each  reach 
are  specified. 

7.  A  complete  mixing  is  assumed  at  all  points  where  a  tributary 
or  wastewater  effluent  enters  a  stream. 

8.  The  design  flow  used  for  determining  the  capacity  of  the 
stream  to  assimilate  wastes  is  usually  the  minimum  average 
consecutive  seven-day  flow  expected  once  in  10  years  on  the 
average  [19]. 

3.3.2   Constraints 

It  is  necessary  to  compute  the  BOD  and  DO  concentrations  at 
the  beginning  and  end  of  each  reach.  In  the  inventory  equations 
that  follow,   the  subscript  r  denotes  a  particular  reach. 

In  general  the  total  flow,  QS  ,  in  reach  r  is  the  sum  of 
the  flow  in  the  previous  reach,  QS   ;  the  tributary  flow  entering 
the  reach,  QT  ;  and  the  wastewater  flow  discharged  into  the 
reach,  QW  . 

Q3r  =  Q3    +  QT^  +  QWr  (10) 

complete  mixing,  the  BOD  concentration  at  the  beginning 
of  each  reach,  EB  ,  is  equal  to  the  sum  of  the  BOD  concentrations 
at  the  end  of  the  previous  reach,  BE   ;  in  the  tributary,  BT  ; 
and  in  the  wastewater  effluent,  BW  ,  times  their  respectivs  flows 
divided  by  the  total  flow. 


Sinilarly,  the  DO  concentration  at  the  beginning  of  each 
reach,  C3r>  can  be  determined  from  the  concentration  at  the  end 
of  the  previous  reach,  CE^j  the  tributary  concentration  CT  j 
and  DO  concentration  in  the  wastewater  effluent,  CW  . 

CEr_l  ^.i  +  CT   QT  +  CW„  qv 

CBr  = __ E i__£  {12) 

The  DO  deficit  at  the  beginning  of  each  reach,  D3  ,  is  the 
difference  between  the  saturation  concentration  CS  ,  and  the 
initial  DO  concentration,  C3  . 

DBr  =  C3r  -  CH  (yj 

The  BOD  concentration  and  CO  deficit  at  the  end  of  each 
reach  can  be  computed  from  the  initial  BOD  concentration  and  DO 
deficit,  using  either  equations  (3)  and  (k)    or  equations  (6)  and 
(7).   The  time,  t,  In  these  equations  is  understood  to  be  equal 
to  the  time  required  for  water  to  flow  from  the  beginning  of 
each  reach  to  the  end  of  that  reach,  T  .   Thus,  BE  ,  the  BOD 
concentration  at  the  end  of  reach  r  and  DE  ,  the  DO  deficit  at 
the  end  of  reach  r,  can  be  determined. 

Using  these  equations  It  la  possible  to  write  constraints 
that  define  the  minimum  allowable  DO  concentration  within  each 
reach . 

The  do  deficit,  Drt,  at  various  points  t  along  each  reach  r 


;ould  be  constrained  to 


.  to  the  maximum 


allowable  deficit  In  that  reach,  Dir'ax. 

Drt  <  ^X   for  various   t  ,   0  <  t  <  Tr  (W 

3y  reducing  the  interval  between  successive  fa,    the  possibility 
of  violation  between  these  points  is  reduced.   If  the  time  of 
flow,  T^,  is  less  than  the  critical  time,  t  ,  only  one  quality 
constraint  is  necessary  for  that  reach,  namely 


0 


(15) 


A  few  trial  solutions  will  enable  one  to  eventually  place 
these  quality  constraints,  equation  (1*0,  at  the  points  having 
the  lowest  DO  concentration.   With  these  constraints,  the  solution 
will  yield  the  maximum  amounts  of  30D  that  can  be  released 
into  each  reach  without  violating  the  standard  for  any  reach. 

Sometimes  an  additional  constraint  may  be  desired  if  each 
treatment  facility  is  required  to  remove  the  same  fraction  of  BOD 
from  the  wastewater  influent.   Such  a  requirement  can  be 
expressed  by  equating  for  each  reach  the  ratios  of  the  BOD 
concentration  released  into  a  reach,  B'.^,  over  the  total  amount 
available,  BWmax. 


BW„      BW" 


r  =   Reach  number 
O  =  Wastewater  treatment  facility 


Fig.  J.      Hypothetical  River  Basin  [19] 


P   =  percentage  treatment  at  r   plant 

Also,  In  general,  F  is  constrained  between  two  specified 
limits. 

i\fn  <   Pr  <  Pjf*  (18) 

3.3-3   Objective  Function 

Objective  function  should  enable  one  to  specify  amount  of 

wastewater  treatment  required  to  meet  at  minimum  cost  a  set  of 

quality  standards  for  the  river  basin.   The  objective  is  to 

minimize  the  total  cost  of  wastewater  treatment,   mathematically, 

it  may  be  stated  as 

N 
MINIMIZE  '   E   C  (P  )  (19) 

r=l  r  r 

where 

C  (P  )  =   the  function  representing  the  total  cost  of 

providing  the  treatment  P  at  the  rth  discharge, 

given  as  an  annual  cost  including  anorlzation 

and  operating  costs. 

The  objective  is  constrained  by  the  quality  standards  as  expressed 
by  (I*)-)  and  (15)  and  the  Inventory  equations  discussed  earlier. 

3.^    DESCRIPTION  OP  SYSTEM 

A  hypothetical  river  basin  shown  In  Pig.  3  was  used  by  Loucks, 
Hevelle  and  Lynn  [19]  to  establish  and  evaluate  various  water 
quality  control  policies.   The  same  basin  with  all  the  necessary 


data  is  used  here.   As  a  matter  of  fact  this  hypothetical  system 
was  derived  from  a  highly  simplified  representation  of  Willamette 
River  in  Oregon.   The  main  contributors  of  pollutents  on  this 
stream  are  municipalities  and  pulp  and  paper  Industries. 

In  this  basin  the  quality  of  water  in  seven  reaches  is 
affected  by  the  amount  of  30D  released  from  six  wastewater 
treatment  facilities.  .It  is  assumed  that  these  facilities  exist 
and  are  currently  removing  a  sufficient  amount  of  BOD  to  satlfy 
the  stream  quality  standards.   It  is  anticipated,  however,  that 
by  1980  the  BOD  load  will  increase  considerably.   This  will 
obviously  necessitate  additional  treatment.   It  is  required  to 
determine  the  additional  treatment  necessary  and  the  minimum  cost 
required  to  meet  the  standards  in  1980. 

Table  1  provides  all  the  necessary  stream  and  wastewater 
data.   Table  2  gives  wastewater  treatment  data.   For  all  treatment 
facilities,  a  minimum  removal  of  35,S  has  been  imposed.   The 
intent  of  this  constraint  is  to  require  each  plant  to  provide  at 
least  primary  treatment,  thus  ensuring  the  abscence  of  floating 
solids  in  the  stream.   The  presence   of  these  solids  is  usually 
considered  objectionable  even  though  they  may  not  reduce  the 
oxygen  concentration  below  the  minimum  acceptable  level. 
.Because  of  technological  difficulties  in  construction  of  facilities 
that  can  remove  over  90 4  of  the  BOD  with  certainty,  $0%   will  be 
assumed  to  be  the  maximum  treatment  required. 

The  treatment  plant  costs  are  usually  convex  within  the 
range  from  J$%   to  90;!  BOD  removal.   A  typical  cost  curve  is  shown 
in  Fig.  4A.   Loucks,  Revelle  and  Lynn  [19]  assumed  that  oqsta  are 
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linear  within  each  segment.  In  this  work,  a  quadratic  fit  was 
obtained  for  treatment  costs  at  each  plant  using  the  standard 
regression  analysis  procedure  [?].  A  typical  quadratic  fit  is 
shown  in  Fig.  43.   The  general  cost  equation  can  be  written  as 

a  +  bP  +  c?2  (20) 

where  F  =  percent  treatment  provided  at  a  plant 
The  coefficients  a,  b  and  c  for  various  plants  are  given  in  Table 
3.  This  particular  aspect  is  a  significant  change  from  [l°]  and 
represents  a  riore  realistic  situation.  The  problem  no  longer  can 
be  solved  by  linear  programming.  It  becomes  one  of  nonlinear 
programming  where  a  nonlinear  objective  function  is  subjected  to 
linear  constraints. 

The  objective  function  can  now  be  written  for  the  river 
basin  illustrated  in  Fig.  3. 

MINIMIZE       Za+bP+cP2  (?1) 

r=1     r         r  r         r  r 

v-15 

The  constraints  must  bound  all  Pr*s,  define  the  J 

final  BOD  and  DO  concentrations  in  each  reach  and  limit  the  ROD 

concentrations  at  the  beginning  of  each  reach  so  that  the  stream 

standard  is  not  violated. 

3.5   SOLUTION 

The  problem  formulated  above  was  a  nonlinear  programming 
problem  with  linear  constraints.   Application  of  the  gradient 
projection  method  seems  quite  justifiable  to  obtain  the  solution 
of  this  problem. 


Two  different  runs  were  made  with  the  data  of  Tables  1  and 
2.   For  each  run,  both  Streeter-Phelps  and  Camp-Dobbins 
formulations  were  used.   Run  1  was  used  to  find  the  optimum 
(minimum)  cost  configuration  of  plants  that  will  just  meet  the 
DO  standards  specified  for  the  stream.   Run  2  duplicated  the 
conditions  of  Hun  1,  except  that  the  minimum  DO  standards  in  each 
reach  had  been  reduced  by  0.5  rag/1. 

A  brief  description  of  different  computational  aspects 
involved  in  solution  by  the  gradient  projection  method  is  not  out 
of  place. 

The  various  tolerances  and  limits  are  to  be  judiciously 
selected.   A  very  small  value  of  £  may  require  a  considerably 
more  computation  time  without  actually  contributing  much  to  the 
improvement  of  functional  value.   In  general,  the  following 
approach  was  used.   Firstly,  using  a  relatively  large  value  of 
£^0.005),  the  solutions  were  obtained  starting  with  different 
initial  points.   Using  this  information,  a  new  starting  point- 
very  close  to  optimal  values  of  control  variables  was  established. 
The  value  of  ^  was  then  considerably  reduced  (.0001)  and  an 
accurate  optimal  solution  was  obtained. 

The  following  values  were  assumed 

£,  =  .005  and  .00  01    MANU  -.=  20       Y    =  0 
£0  =  -°05  MXHN  =   3       ,,i    =  0 


The  problem  being  one  of  minimiiiatlor 


gradient  projection  maximization  algorithm  by  maximizing  the 
negative  of  the  original  objective  function.   The  constant  U 
appearing  in  the  objective  function  given  by  (21)  was  not 
considered  in  optimization.   Thus, 


while,  the  actual  total  minimum  cost  is  given  by 

7 

F'  =   2  a   -  maximum  (-F)  (23) 

r=l  r 

Five  different  initial  points  were  used  for  each  case. 
This  information  is  shown  in  Tables  U-   through  7  along  with  details 
about  execution  times,  number  of  iterations  and  number  of 
functional  evaluations.   The  values  of  lip  gjl  at  optimal  are  also 
given.   Tables  8  through  11  indicate  the  convergence  rates 
obtained  for  certain  specific  cases.   The  same  information  is 
illustrated  in  Fig.  5  for  one  selected  case. 

The  optimal  (minimum)  cost  solutions  obtained  for  each  case 
are  presented  in  Tables  12  through  15.   These  tables  also 
provide  the  information  regarding  the  maximum  amount  of  EOD  that 
can  be  released  to  the  stream  and  the  resulting  CO  concentrations. 
Tables  16  and  17  indicate  the  linear  programming  solutions 
obtained  by  Loucks,  Revelle  and.  Lynn  [19]. 

3.6   DISCUSSION 

The  model  presented  in  this  work  can  be  used  to  determine 
the  minimum  total  cost  associated  with  any  particular  set  of 
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%   BOD     Annual   Maximum  Minimum  Minimum 

BOD  Release  DO  in  minwMo 

Removal   Cost  (»   to  Stream  Reach  DO  in  Reach 

(mg/l/day)  (mg/i)     (mg/1) 


1 

66.5 

0 

S3. 00 

9.50 

7.00 

2 

63.0 

631,200 

152.70 

7.61 

7.50 

3 

48.0 

171,921 

124.80 

8.15 

7.00 

'+ 

90.0 

695.350 

lWj-.OO 

6.01 

6.00 

5* 

- 

- 

- 

7-95 

6.50 

6 

90.0 

890,057 

218.00 

6.91 

6.00 

7 

63.8 

873,718 

101.00 

4.0  0 

4.00 

Total      3,262,246     823. 50 
3  wastewater  effluent  is  released  into  Reach  5. 


l^Tm/V'   3°1Uti0n  When  D0  standards  Are  Reduced 
Streeter-Phelps  Formulation 


Reach  No. 

%   BOD 
Removal 

Annual 
Cost  (|) 

Maximum 
BOD  Release 

to  Stream 
(mg/l/day) 

Minimum 
DO  in 
Reach 
(ns/D 

Minimum 
Allowable 
CO  in  Reach 
(mg/l) 

1 

65 

0 

85.90 

9.50 

6.50 

2 

56 

575,890 

178.20 

7.41 

7.00 

3  . 

^3 

16^,365 

136.20 

8.  04 

6.50 

4 

90 

695,360 

144.00 

6.  01 

7.  86 

5.50 
6.00 

6 

90 

890, 042 

218.00 

6.  80 

5.50 

7 

52 

748,963 

132.60 

3.  50 

3.50 

Total      3,074,620     894.90 
Mo  wastewater  effluent  is  released  into  Reach  5. 


Reach  No. 

,?  BOD 
Removal 

Cost  (3) 

BOD  Relase 
to  Stream 
(mg/l/day) 

I-'iniaun 
DO  In 
Reach 
(mg/1) 

Minimum 
Allowable 

DO  in  Reach 

(Dlg/1) 

1 

66 

0 

34.0 

9.50 

7.00 

2 

60 

607,500 

162.4 

7.66 

7.50 

3 

W.4 

169,086 

128.5 

8.32 

7.00 

4 
5* 

90 

695,368 

144.0 

6.04 
8.19 

6.00 
6  .50 

6 

90 

890,020 

218.0 

7.17 

6  .00 

7 

62.8 

860,762 

103.8 

4.12 

4.00 

Total       3,222,736    840.7 
)  wastewater  effluent  is  released  into  Reach  5. 


r'lninura  Cost  Solution  When  DO  Standards  Are  Reduced 
by  0.5  mg/1 
Camp-Dobbins  Formulation 


BOD     Annual    Maximum  Klnlmum  Kinimua 

-,  '  „  .    ,,,      SOD  Release  DO  in  Allowable 

Moval  Cost  ($)   to  Stream  Reach  DO  in  Beach 

(rag/1/day)  (r.g/1)     (mg/1) 


0 

86.92 

9.50 

6.50 

560,71^ 

I87. 95 

7.^7 

7.00 

162,607 

139.82 

8.22 

6.50 

695,368 

1W.00 

6.0^ 

5.50 

- 

- 

8.10 

6.00 

890,022 

218.00 

7.06 

5.50 

730,226 

138.67 

3.61 

3.50 

Total      3,038,937     915.56 
3  wastewater  effluent  is  released  into  Reach  5. 


Table  16.   Minimum  Cost  Solutlo.. 

Oxygen  Standards  by  Loucks  et.al.  ^19] 


intaininR:  Dissolved 


Reach  No  *   BODReleased  Minimum  Minimum 

Reach  No.  to  Stream     DO  in  Allowable 

Removal  Cost  (|)   (mg/l/day)     Reach  DO  (mg/1) 
(mg/1) 


608,500 
170,000 

690,000 

900,000 
902,000 


Total       3.270,500       8^9 
5  wastewater  effluent  is  released  into  Reach  5. 


Table  17.   Minimum  Cost  £ 
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Keach  No.  to  stream     20  in    Allowable 

Removal   Cost  (*)   (mg/l/day)     Reach    DO  SI) 


522,000       20^        7.3  7.0 

170,000       120        8.2  6.5 

690,000      1W       6.0  5.5 

7.9  6.0 

900,000       218        7.0  5.5 

690,000       139        3.c  ,  « 
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minimum  allowable  DO  concentrations  in  a  river  basin.   This 
model  also  provides  an  useful  tool  to  determine  sensitivity  of 
both  the  cost  and  the  actual  minimum'  DO  concentrations  in  reaches 
to  changes  in  the  minimum  allowable  DO  concentration  In  any 
particular  reach. 

As  the  data  for  the  problem  was  taken  from  the  work  of  Loucks, 
Hevelle  and  Lynn  [19].  a  comparison  of  the  results  is  not  out  of 
place.   Discussion  here  is  confined  to  only  Gamp-Dobbins 
formulation,  but  the  same  applies  to  Streeter-Phelps  formulation. 
The  results  of  Loucks,  Hevelle  and  Lynn  are  summarized  in  Tables 
16  and  1?.   These  correspond  to  Tables  Ht  and  15,  respectively 
of  this  work.  , 

Note  that  whereas  Plants  2  and  3  provided  60%   and  b€>.5% 
treatments  (Table  1'+),  these  plants  provided  55;s  and  50;J  (Table 
16)  in  the  linear  programming  solution  presented  by  them.   But  it 
is  worth  noting  that  the  sum  of  costs  at  Plants  2  and  3  was 
nearly  identical  in  both  cases.   Note,  also,  that  the  costs  for 
the  two  solutions  were  nearly  identical,  this  solution  costing 
:;',h7 ,7k'ir   less  than  the  linear  programming  solution.   This  small 
difference,  may  well  stem  from  the  fact  that  they  used  linearized 
cost  curves  while  a  quadratic  fit  was  used  in  this  solution. 

Table  1J  presents  the  optimal  (minimum)  cost  solution  when 
the  minimum  allowable  DO  concentration  is  reduced  by  0.5  rag/1  in 
each  reach.  This  reduction  in  DO  standards  results  in  an  annual 
cost  savings  of  5183,800  or  about  6  percent  of  the  total  cost. 

From  Tables  14  and  15 ,  it  is  clear  that  the  DO  standards 
for  Reaches  2  and  7  dictate  the  actual  concentrations  in  the 


entire  basin.   In  other  words,  a  reduction  in  the  minimum  DO 
concentrations  in  any  but  Reaches  2  and  7  would  neither  decrease 
the  minimum  total  cost  nor  the  actual  DO  concentrations.   Table 
15  shows  that  once  the  DO  standards  have  been  reduced  by  0.5  rag/1 
in  each  reach,  only  the  Reach  ?,  can  be  regarded  as  critical 
and  hence  it  can  be  said  to  determine  the  required  treatment, 
and  therefore,  the  cost,  throughout  the  basin. 

It  is  evident  from  the  results  that  a  change  in  the  minimum 
CO  concentration  standards  in  several  reaches  may  have  no  effect 
on  the  DO  concentrations  in  these  reaches.  Conversely,  a  change 
in  the  minimum  allowable  DO  concentration  in  a  single  reach  may 
affect  the  DO  concentrations  in  every  other  reach. 


A  decrease  in  the  minimum  DO  c 


ty   0.5  mg/1  in 


each  reach  of  the  hypothetical  basin  only  reduces  the  actual 
minimum  DO  concentrations  by  0.2  mg/1  in  Reach   2,  0.1  mg/1  in 
Reaches  3,  5  and  6  and  0.5  mg/1  in  Reach  7.   The  change  in  annual 
benefits  resulting  from  these  lower  oxygen  concentrations  can  be 
compared  to  the  annual  cost  savings  of  3185,875  in  order  to 
determine  the  desirability  of  this  reduced  standard. 

In  many  basins,  under  present  legislation,  uniform  standards 
in  terms  of  a  final  percentage  removal  of  waste  material  are 
imposed  upon  polluters.   If  such  a  standard  is  applied  to  the 
every  facility  in  this  basin,  it  is  evident  that  every  facility 
will  have,  to  remove  the  same  percentage  as  is  required  by  the 
treatment  facility  on  Reach  V  namely  90*.   Anything  less  than 
90,1  waste  removal  from  the  effluent  entering  Reach  h   would  result 
in  a  minimum  DO  concentration  less  than  the  minimum  allowable  DO. 


The  same  conclusion  can  not  be  drawn  from  noting  that  90,£ 
treatment  is  required  on  Reach  6.   Since  in  this  case  there 
exists  some  control  over  the  concentration  of  waste  and  DO  in 
the  water  entering  that  reach. 

The  gradient  projection  method  did  not  encounter  any  trouble 
in  solution  of  the  problem.  The  efficiency  of  method,  of  course, 
depends  upon  a  judicious  choice  of  various  limits  and  tolerances. 


'   :  AiiD  CONTROL  0?  RbS^RVOIR  SYSTEMS 
4.1   INTRODUCTION 

The  aathematioal  models  which  are  used  to  describe 

resource  systems  often  contain  nonlinear  mathematical  relationships 
that  are  difficult  to  analyze  and  optimize.   Furthermore ,  a  large 
water-resource  system  is  generally  a  multidimensional  problem. 
Pioneering  work  in  water-resource  systems  and  optimization 
analysis  has  been  carried  out  in  the  Harvard  Water  Program  [l6, 
20].   locations  of  recent  research  [16]  imply  that  simulation 
is  still  being  used  in  the  del  died,  final-stage  optimization  of 
a  given  water-resource  system. 

Hall  and  others  [13,  lfc]  Were  the  flrBt  to  propOGe  the 
application  of  dynamic  programming  to  the  optimization  of 
reservoir  systems.   Meyer  [22]  successfully  used  dynamic 
programming  approach  in  optimizing  the  operation  of  a  multiple- 
purpose  reservoir.   Dynamic  programming  has  the  advantage  of 
effecting  the  decomposition  of  a  highly  complex  problem  into  a 
series  of  far  less  complex  problems.   However,  due  to  the 
dimensionality  difficulties,  this  procedure  can  not  be  extended 
to  more  practical  problems. 

Linear  programming  is  not  generally  useful,  because 
individual  reservoir  utility  functions  are  in  general  definitely 
nonlinear.   Quite  often,  many  of  the  restrictions  on  the  operation 
of  a  water-resource  system  are  linear.   This  situation  arising 
from  a  nonlinear  objective  or  utility  function  subjected  to 


linear  constraints  provides  an  useful  field  for  nonlinear 
programming  techniques  suoa  us  the  gradient  projection  method. 

Specifically,  the  gradlSnt  projection  method  discussed 
'-   :'   ',  Is  applied  here  to  solve  two  simple  problems.   The 
stress  is  on  the  method  that  ho;,  a  general  nonlinear  objective 
can  be  treated  without  going  into  tedious  details  of  H,   , 
Both  hypothetical  systems  are  taken  from  [20]  along  with  all  the 
relavent  data  needed  to  solve  them. 

Initially,  the  artificial  system  from  which  these  two 
hypothetical  sytems  are  derived,  is  described.   This  includes  a 
brief  description  of  streamflow  data  a,d  certain  basic  assumptions 
regarding  irrigation  and  water  power.  Following  this,  both 
models  are  described  and  solved  by  the  gradient  projection 
method . 

4.2    DESCRIPTION  Op  SYSTEM 

The  simplified  system  ba.od  on  the  Clearwater  Silver  Basin 
in  Idaho  is  described  in  detail  in  the  Harvard  Water  Program 
[20].   For  this  system,  at  least  one  of  each  major  Kind  of  output 
of  a  water-resource  system  was  developed.   These  were,  a  with- 
drawal -  consumptive  use;  a  nonwlthdrawal ,  essentially  non- 
consumptive  use;  and  a  retardation  or  Withholding  use.   For  these 
purposes,  the  irrigation  of  crops,  the  development  of  water 
power  and  reduction  of  flood  damage  were  considered. 
4.2.1  Physical  Layout  of  the  System 

The  physical  layout  of  the  system  chosen  for  development  is 


Point  gL  db     Power  plant 

Power  plant  G  □  A 

§,      Reservoir 
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Sketch   of  the  Simplified   River-Basin  Sy8t« 


shown  in  Fig.  6.  There  are  fo  ,  |  ,  <au  lding  reservoirs  A,  B,  and 
C  and  Di  two  power  plants,  one   at  reservoir  B  and  other  at  point 
G.   Also  there  is  an  irrigation  diversion  dan  at  point  E.   Only 
reservoirs  A  and  B  can  provide  releases  for  irrigation. 
Reservoir  D  lies  on  a  tributary  stream  which  joins  the  main  stem 
at  point  P.   Irrigable  areas  lie  on  both  banks  of  the  main  river 
downstream  from  reservoirs  A  and  B.   Return  flow  reenters  the 
river  channel  below  E  but  upstream  from  reservoir  C;  none  of 
it  can  be  reused  within  the  irrigated  areas.   A  flood  damage 
zone  is  situated  just  below  point  G. 

The  additional  details  have  been  purposely  avoided  in  the 
development  of  the  system  because  these  may  unnecessarily 
complicate  operation  studies,  the  computer  programming  and  the 
general  analysis. 

k.Z.Z       Streanflow  Data 

A  detailed  description  of  the  str  mflow  data  is  gin  „  in 
[20].  Table  13  identifies  the  magnitudes  of  mean  monthly  flows 
observed  at  point  E.  The  observed  pattern  is  typical  of  the 
hydrology  of  a  catchment  area  in  which  the  melting  of  winter 
snows  produces  large  spring  runoffs,  followed  by  low  summer  and 
fall  discharges.  Such  basins  are  found  in  wide  regions  of  the 
wastern  United  States.   Data  for  other  points  also  indicate  the 
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Lve  use  and  diversion  requirement,  for  irrigatio 
itological  data  and  irrigation  practices  in  the 


region.   The  region  is  semiarid  with  a  moderately  long  growing 
season  of  2G.5  days. 

The  unit  irrigation-diversion  requirement  was  assumed  to 
be  5.0  ft(op  5.0  acre  ft  per  acre)  yearly.   This  t*  distributed 
by  months  as  shown  in  Table  19. 

The  estimate  of  return  flow  from  irrigation  diversion  was 
based  on  certain  assumptions.   Mo  water  would  be  lost  by 
evaporation  from  drainage  and  wasteway  channels,  that  no  return 
flow  would  be  consumed  on  nonlrrlgated  land,  and  that  no  water 
escape  into  neighbouring  basing.   Return  flow  would  therefore 
equal  the  difference  between  the  irrigation  diversion  requirement 
and  the  set  consumptive  use  of  irrigation  water.   In  general, 
it  comes  to  about  50  to  60  percent  of  total  irrigation  diversion 
requirement.   The  monthly  distribution  of  this  return  flow  in 
years  of  full  irrigation  supply,  shown  in  Table  20,  is  in 
accordance  with  observations  in  several  Bureau  of  Reclamation 
projects. 

Irrigation  output  of  the  system  was  set  at  6  x  106  acre  ft. 
Unit  annual  gross  irrigation  benefits  were  taken  as  decreasing 
from  about  %   6.50  per  acre  feet  at  low  levels  of  development 
to  about  %   5.00  per  acre  feet  at  maximum  development.  Figure  7 
shows  the  unit  gross  benefit  function. 

The  losses  cccuring  from  shortages  are  not  discussed 
bere.   The  estimates  of  capital  costs  and  costs  of  operation, 
maintainance.  and  replacement  (OMH)  for  irrigation-diversion 
works  are  given  in  [20].   They  are  based  on  data  from  Bureau 
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«ont. 

~  Flow 
(10Z  acre  ft) 

Percentage  of 
Mean  Annual  Runoff 

January 

1,698 

3.00 

February 

1.778 

3.10 

March 

3,066 

5.40 

April 

3,939 

15.80 

Hay 

18,100 

32.00 

June 

12,618 

22.10 

July 

3,469 

6.20 

August 

1,079 

1.80 

September 

845 

1.50 

October 

1,264 

2.20 

November 

1,813 

3.30 

December 

1,989 

3.60 

Total 

56,658 

100.00 

Assumed  Monthly  Distribution  of  Annual 
Irrigation  Diversion  Requirement  [20] 
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Fig.  7.   Assumed  Unit  Gross 


Irrigation  Benefit  Function  [20] 
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Assume-!  Monthly  Distribution  of  Annual 
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Supply  [20 J 
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Table  21.   Assumed  Monthly  Distribution  of  Annual 
Energy  Requirement  [20] 


November 
December 
January 
February 
March 


h.9  -  y  -  : 
\l.8  +  Y  +  ; 

U-9    -   Y    -   Z   -    0.2751 

\3.9  +  X  +  Z  -  0.1251 


[20] 


i   of   Configuration  Use'-i    in  Two-I 


of  Reclamation  projects. 
'1.2.'!   Water  Power 

The  power-requirement  of  the  system  was  selected  to  be 
representative  of  a  diversified  agricultural-industrial,  rural- 
urban  economy.   The  assumed  monthly  distribution  of  the  load 
is  listed  in  Table  21.   The  maximum  yearly  target  output  for 
energy  for  the  system  was  set  at  k   x  109  kw  hr. 

Unit  energy  benefits  were  assumed  to  be  constant  for  all 
levels  of  putput,  namely.  9  mills  per  kw  hr.  of  firm  energy 
and  1,5  mills  per  kw  hr.  of  nonfirm  energy.   Assumed  capital 
costs  an4  0MB  costs  for  power  plants  are  given  in  [20]. 

*.3   A  MODEL  WITH  TWO  SEASONS  AND  PREDICTABLE  HYDROLOGY 

This  example  deals  with  a  relatively  simple  configuration 
of  uses  and  installations  under  conditions  in  which  the  pattern 
of  inflows  repeats  itself  each  year  with  certainty  so  that 
overyear  storage  is  not  required.   This  is  the  simplest  problem 
that  could  be  devised  while  retaining  enough  suvstance  to 
present  some  challenge. 

The  configuration  shown  in  Fig.  8  has  been  abstracted  from 
simplified  river  basin  discussed  earlier.   This  configuration 
retains  two  reservoir  sites,  an  irrigation  project,  and  a  run- 
of-the  river  power  plant. 

Most  of  the  data  are  contained  in  the  figure.   The  numbers 


70 
in  braces  give  mean  flows  In  the  river  in  two  seasons  of  the 
year  the  wet  season  (top  figure)  and  the  dry  (bottom  figure). 

The  first  topographic  feature  the  river  encounters  is  the 
reservoir  B,  its  active  capacity,  denoted  by  X,  is  one  of  the 
unknowns  of  the  problem.   The  amount  of  water,  Y,  will  be 
retained  in  the  wet  season  and  released  in  the  dry. 

Just  below  the  confluence  of  the  west  branch  with  the  main 
stem  an  irrigation  diversion  canal  takes  off  water  to  an  Irrigated 
area  lying  to  the  east'.'  The  total  amount  of  irrigation,  I,  is 
the  second  unknown  of  the  problem.   It  is  assumed  that,  whatever 
I  may  be,  k2.5   percent  of  irrigation  water  must  be  provided  in 
the  wet  season  and  57.5  percent  in  the  dry.   The  resulting  return 
flows  from  the  irrigated  area  are  assumed  to  be  15  percent  of  I 
in  the  wet  season  and  ^5  percent  in  the  dry. 

It  is  also  necessary  to  determine  the  usable  capacity  of 
reservoir  C,  denoted  by  Z.   The  fourth  and  final  variable  is 
the  energy  output  of  power  plant,  E.   Half  the  annual  output  of 
energy  generated  is  assumed  to  be  required  in  the  wet  season 
and  half  in  the  dry. 

h.J.X       Constraints 

The  first  group  of  constraints  requires  simply  that  none 
of  the  four  decision  variables  be  negative. 


The  second  group  of  constraints  states  that  the  flows  in 
all  reaches  of  the  system  must  be  nonnegative.  From  the  map 
these  constraints  can  be  written  down. 

3-3-2  >  0 

3.9  -  X     -  0.4251  >  0 

1.8  +  X     -  0.5751  >  0 

3.9  -  Y  -  Z  -  0.2751  >  0 

Although  there  are  other  six  flow  constraints  involving  decisio: 
variables,  these  are  autometically  satisfied  if  above  four 
constraints  are  satisfied. 

The  third  group  of  constraints  asserts  that  the  flow  at  the 
power  plant  must  be  adequate  in  both  the  wet  and  the  dry  seasons 
to  generate  the  amount  of  power  that  has  been  decided  on. 

The  technical  relationship  between  flow  and  energy  output 
is  taken  to  be 

E  =  O.l^P 

E  *   energy  generated  in  any  period  in  109  kw  hrs. 
and    P  =   flow  through  turbines  in  106  acre  feet 

As  equal  amount  of  energy  is  reguired  in  both  n-aasons,  the  two 
power  constraints  are 

6.9  -  Y  -  Z  -  0.2751  >  C.5E/0.133  =  3-^7E  V 

and 

3.9  •[■  Y  -  Z  -  0.1251  >  0.5E/0.144  =  3.47E  (1, 


t  i  . ,  these  power  oomstralnts  for  wet  i 
seasons,  respectively  oan  be  written  as 

Y  +  Z  +  o.2?5I  +  J.h'ni  <  6.9  (H) 

-Y  -  Z  +  0.1251  4  3.^72  <  3-9  (12) 

The  design  sought  is  assumed  to  be  the  one  which,  while 
satisfying  these  constraints,  yields  the  greatest  possible  present 
value  of  net  benefits. 

'+•3.2   Objective  Function 

The  objective  function  to  be  maximized  is  [20] 

tt  ■•  31(E)  +  32(I)  -  Ca(Y)  -  C2(Z)  -  C3(E)  -  c4(I)     (13) 


n     =  the  present  value  of  net  benefits  in  106  dollars 
B1(E)  =  the  present  value  of  an  output  E  x  109  kv;  hr  per 

year  in  10   dollars 
E2(I)  =  the  present  value  of  an  irrigation  supply  of 

I  x  10   acre  ft  per  year  in  106  dollars 
C^Y)  -,  the  capital  cost  of  building  reservoir  3  to  capacity 

Y  in  106  dollars 
C2(Z)  =  the  capital  cost  of  building  reservoir  C  to 

capacity  Z  in  106  dollars 
C3(E)  =  the  capital  cost  of  build  in?  the  power  plant  to 

capacity  E  per  year  in  106  dollars 

C4(I)  =  the  capital  cost  of  building  the  irrigation  system 


to  capacity  I  per  year  in  10  dollars 

The  data  for  all  these  functions  are  given.   Firstly,  capital 
cost  functions  arc  discussed. 

c (r)  .  4jy/(i  +  0.2Y)  (14) 

cg(z)  =  4?z/(i  +  o.3z)  (15) 

C  (E)  =  20. 6E  -  E2  (16) 

It  can  be  seen  that  the  larger  the  reservoir  capacity,  the 

higher  is  the  cost.   E  is  obviously  restricted  to  a  higher  value 

pf  4  x  105  kw  hr  as  indicated  earlier.   The  derivation  of  C. (I) 

is  a  bit  complicated  because  of  the  assumption  that  only  3  x  10 

can  be  taken  for  irrigation  without  pumping.   A  pumping  plant  is 

required  if  more  than  this  amount  is  required  for  irrigation. 

The  data  to  be  assumed  are  as  follows. 

The  basic  cost  of  the  diversion  works  is  $4,500,000  plus 

$44,000,000  per  10  acre  ft  of  Irrigation  water.   For  I  > 

6 
3  x  10  acre  ft,  a  pumping  olant  must  be  constructed  at  a  cost 

6 
of  $500,000  plus  $20,000,000  per  10  acre  ft  of  water  to  be 

pumped.  Now  if,  I  denotes  non  pumped  portion  and  I  denotes  the 

pumped  portion,  then  the  capital-cost  function  for  the  irrigation 

works  becomes  - 

C  (I)  o  441  +  641  +  4.51*  +  0.51*  (17) 


Ix<3 


I  =  J    if   I  >  0 

4   -  1    if   J2  >  ° 
The  next  step  is  to  formulate  the  benefit  functions  3  (E) 
and  Bg(I).   The  calculation  of  each  of  these  proceeds  in  four 


Express  annual  gross  benefits  as  a  function  of  E  or  I . 

Express  annual  operation,  malntalnance,  and  replacement 

(OKH)  costs  as  a  function  of  E  or  I. 

Compute  annual  net  benefits  by  subs  traction. 

Compute  the  present  value  of  net  benefits  by  applying 

an  appropriate  present-value  factor. 

Assuming  a  planning  period  of  50  years  and  a  discount  rate 
of  2  1/2  percent,  a  present  value  factor  of  28.^  is  obtained. 
That  is,  under  these  assumptions  the  present  value  of  a  net 
benefit  of  $1  per  year  for  50  years  is  $28 .40. 

B  (E)  is  derived  as  follows.   Assuming  that  all  the  energy 
Is  on  demand,  a  price  of  9  mills  per  kw  hr  can  be  assumed. 
Therefore,  the  gross  benefits  in  10  dollars  are  9E.   OKR  costs 
are  assumed  to  be  0.2E  per  year.   This  leads  to  annual  net 
benefits  of  8.8E  from  electric  power.   It  follows  that  the 
present  value  of  electric  power  operations  is 

B1(E)  =  28.^  x  8.8S  =  250E  (18) 

The  calculation  of  B  (I)  is  comparatively  more  complicated 
because  of  following  reasons.   The  introduction  of  pumping  plant 
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causes  discontinuity  and  also  marginal  value  of  irrigation  water 
can  not  be  regarded  as  constant.  The  datum  assumed  for  computing 
the  gross  benefit  from  irrigation  is 

Marginal  gross  benefit  =  2.1  +  3.2/(1  +  0.21)       (19) 

Integrating  this  from  0  to  I,  total  gross  benefits  in  10  dollars 
can  be  obtained. 

Total  gross  benefitsi   2.11  +  36.8log(l  +  0.21)        (20) 

The  OMR  costs  can  be  assumed  to  be  0.51.  +  l-56l2.   Thus, 

Annual  net  benefits  =  1.6l,  +  0.541  +  36.8  log  (1  +  0.21) 

(21) 
Finally,  applying  the  present-value  factor  1 

B2(I)  ^  45.41^  +  15.3I2  +  10451og(l  +  0.21)  (22) 

The  objective  function  can  now  be  computed  by  adding  expressions 
(I»0,  (15),  (16).  (17),  (18)  and  (22).   Hencei 

rr  =  229. 4E  +   E2  +  1.4l1  -  48.7I2  +  1045105(1  +  0.21) 

-  ^-51^  -  3.. -".I*  -  43Y/(1  +  0.2Y)  -  4725/(1  +  0.3Z)     (23) 

4.3.3   Solution 

This  problem  requires  finding  the  maximum  of  a  nonlinear, 
and  indeed  discontinuous,  function  of  some  decision  variables 
that  are  related  by  a  number  of  linear  constraints.   Though 
simple,  such  a  formalization  is  appropriate  for  the  initial 
analysis  of   many  water-resource  design  problems. 

Apart  from  the  complexity  of  the  objective  function,  a 


this  sort  can  be  BOJ  forwardly  by  the  well 

.ational  technique  of  linear  programming.   As  the 
objective  function  is  'separable',  the  replacement  of  nonlinear 
expressions  by  linear  segments  is  possible.  This  will,  of  course, 
result  in  introduction  of  many  new  variables.   Such  a  procedure, 
though  useful  for  an  ultimate  application  of  linear  programming 
is  not  recommended  because  that  will  considerably  increase  the 
computation. 

Instead  of  attempting  the  simplification  of  the  objective 
function,  the  applicability  of  the  gradient  projection  method  was 
tested.   A  brief  description  of  computational  aspects  follows. 

The  various  tolerances  and  limits  required  for  the  applica- 
tion of  the  method  were  selected  in  the  same  manner  as  discussed 


:thod   were 

selected    in   the   sami 

s  manner  s 

;ifically, 
1005 

the  following  value; 
MXNU  =   20 

"T 

)05 

KXRN  =      3 

^max 

)05 

f>max  =     1 

Tmax 

The-  problem  being  one  of  maximization  was  solved  by  the 
gradient  projection  maximization  algarithm  directly. 

Number  of  runs  were  tried,  each  having  different  Initial 
values  for  control  variables.   In  selecting  these  initial  stari 
values,  one  has  to  be  quite  reasonable.   If  values  selected  are 
quite  far  off  the  optimal  values,  a  constraint  violation  may 
result.   This  would  require  a  re-inversion  to  continue  and  mor< 
the  number  of  re-inversions,  less  accurate  is  the  solution. 


Therefore,  the  following  a;  opt„d.   ?or  each  variable  a 

range  of  value-  was  selecte",   ftn  Initial  point  then  was  selected 
having  bhe  value  for  each  variable  within  its  specif led  range. 
This  approach  turned  out  to  be  more  succeful  than  one  of  giving 
uniform  values  to  all  variables.   Results  of  these  different  runs 
are  summarized  in  Table  22. 

The  optimal  solution  is  presented  in  Table  2.3.  The  solution 
to  the  same  problem  was  obtained  in  [20],  by  the  method  of  chorda! 
approximation.   Table  Zh   indicates  the  details  of  that  solution. 

The  gradient  projection  method  did  not  encounter  any  diffi- 
culty in  reaching  the  optimal  solution,  though  the  objective 
function  involved  logarithmic  expressions.   The  method  of  chordal 
approximation  reduces  the  accuracy  of  the  overall  no-5  el,  because 
approximate  functions  are  introduced.   Of  course,  the  loss  of 
accuracy  can  be  compensated  for  by  employing  a  finer  grid  for 
the  straightline  approximations.   No  doubt,  the  better  approach 
will  be  to  treat  the  nonlinearity  directly  and  this  is  done 
quite  efficiently  by  the  gradient  projection  method. 

k.''\        A  MODEL  WITH  KOBE  SEASONS  A  P  I     3W  TE  TDIiOLOGY 

The  problem  in  the  previous  section  indicates  that 
mathematical  programming  is  quite  helpful  in  finding  the  optimal 
designs  and  operating  procedures  for  a  fairly  comulex  water- 
resource  system  provided  a    short  segment  of  time  can  be  cor.sidere-" 
In  isolation.   This  limitation  arises  from  the  fact  that  the 
programming  computations  rapidly  become  more  difficult  as  the 
number  of  constraints  to  be  handled  increases.   An  increase  in 
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Fig.  9.   Sketch  of  Configuration  used  in  Four-Period  Problem  [20] 


number  of  periods  increa;  of  constraints. 

In  the  present  section  a  problem  involving  four  periods 
and  requiring-  an  overyear  storage  is  considered.   This  latter 
feature  adds  considerable  complication  to  the  analysis.   Such 
a  situation  is  found  in  regions  in  which  natural  inflow  is  likely 
to  be  deficient  in  some  years,  so  that  overyear  storage  is 
required  for  efficient  operation.   Unpredictability  of  flows  is 
not  considered.   In  other  words,  the  problem  is  of  designing 
a  system  that  makes  efficient  use  of  a  sequence  of  unequal  but 
predictable  flows. 

The  configuration  shown  in  Fig.  9  is  a  modification  of  one 
dealt  with  in  the  previous  section.   In  order  to  emphasize  the 
effect  of  more  periods,  certain  complicating  features  of  last 
section  are  dropped.   The  schematic  representation  indicates  the 
reservoir  C  has  been  suppressed  and  the  irrigation  supply  has 
been  set  at  3  x  106  acre  ft  per  year. 

h.h.l       Constraints 

The  first  set  of  constraints  requires  that  the  volume  of 
water  released  from  the  reservoir  must  be  sufficient  to  meet  the 
period's  irrigation  demand,  where  the  latter  is  a  preassigned 
proportion,  kfc,  of  the  annual  irrigation  demand.   This  requirment 

%  +  f2t  *   V     *  "  !•  ?-    3'  *  (2*> 

where  f^  denotes  the  flow  from  the  tributary  which  joins  the 
main  stem  just  above  irrigation-diversion  canal.   The  natual  flows 
are  indicated  in  the  figure  while  values  of  k,  are  given  in 


Enevsy  demand 


Table  25. 

The  second  constraint  is  that  the  volume  of  water  released 
during  any  period  can  not  kxcpJ  the.  contents  of  the  reservoir 
at  the  beginning  of  the  period  plus  the  flow  in  reservoir  during 
the  period.   Let  S^   denote  the  contents  of  the  reservoir  at  the 
beginning  of  period  I.   Then  the  constraint  is 

at  -  St  *  flt  •       t  '   1'  2-  3.  *  (25) 

where  f.   is  the  preassigned  natural  flow  into  the  reservoir 
during  time  period  t. 

The  third  constraint  is  that  the  contents  of  the  reservoir 
at  the  beginning  of  any  period  can  not  exceed  the  amount  left 
over  from  the  previous  period,  or 

St  *  Vl  +  fl,t-l  -  Vl       t  2,  3.  4  (26) 

Also  it  is  necessary  to  include  a  constraint  which  makes  sure 
that  the  contents  of  the  reservoir  at  the  end  of  any  period,  can 
not  exceed  the  capacity  of  the  reservoir,  or 


l,t 


■■   1.  2,  3.  *» 


The  last  two  requirements  ensure  that  the  contents  of  the 
reservoir  at  the  beginning  of  any  period  do  not  exceed  its 
capacity,  hence  it  is  not  listed  separately. 

The  last  constraint  is  related  to  power  generation.  The 
flew  of  water  past  the  power  plant  must  be  sufficient  to  meet 
the  requirements  of  power  generation.   As  in  the  previous  sect: 


it  is  assumed,  that  9.5  xlo"  acre  ft  of  water  are  required  to 
generate  1  x  10'  kw  hr  of  electric  enorzy.      The  flow  available 
at  the  power  plant  is  the  sun  of  the  flow  past  the  irrigated 
area,  the  return  flow  frps  the  irrigated  area,  k'  I  (where  the 
return  flow  coefficients  are  given  in  Table  25)  and  the  natural 
flow  from  the  eastern  tributary,  f   .   Therefore, 

at  '•  fzt  ~  (kt  "  ki)J   +  f3t  ?-  6-95Ef   *  -  1.  2-  3,  k 

(28) 
where  E^  is  equal  to  a  specified  proportion,  C.  (shown  in  Table 
25)i  of  the  annual  energy  output. 

k.h-.Z        Objective  Function 

Since  the  quantity  of  irrigation  is  prespecif led,  it  no 
longer  enters  the  objective  function.  The  objective  function 
now  consists  of  only  two  fcermsi  the  capital  cost  of  construe  kin:; 
the  reservoir  3  and  the  present  value  of  the  net  hydrodectric 
benefits.  The  expressions  used  in  previous  example  are  used 
here.   Thus 

tt  =  229. U-E   +  E2  -  43Y/U  +  0.2Y)  (29) 

The  purpose  is  to  obtain  the  optimal  value-  of  control 
variables  afct  St>  Y  and  E  which  satisfy  the  constraints  (24) 
through  (27). 

k.h*3       Solution 

Again,  one  approach  to  solve  this  nonlinear  programming 

problem  involving  linear  constraints,  would  be  by  replacing  the 


8? 

nonlinear  expressions  by  linear  segments.   If  a  problem  would 
have  been  linear,  a  method  that  drastically  reduces  the  number  of 
constraints  that  have  to  be  handled  at  any  one  time  can  be  used. 
This  is  the  decomposition  principle  developed  by  Dantzig  and 
Wolfe.  [3].   Such  an  approach  is  very  useful  for  a  large  problem 
and  is  not  required,  here. 

The  gradient  projection  method  was  applied  to  solve  this 
problem.   The  following  values  for  different  limits  and  tolerances 

£x  ■  .0005  MXNU  =20         Y    =  0 

£2  =  .005  MXBN  =3        -«    =  0 

fc3  =  .0005        B^^  =1       y    =  lc 

The  problem  being  one  of  maximization,  the  gradient 
p  'jection  maximization  algorithm  was  applied  directly  without 

A  close  look  at  the  details  of  the  problem  indicates  that 
it  is  not  necessary  to  have  both  St  and  a  in  the  optimization. 
Inequality  (26)  is,  in  fact,  an  equality  and  hence  knowing  one  of 
these,  other  can  be  determined.   Therefore,  in  actual  procedure 
&%   were  not  considered  but  were  derived  from  values  of  S  . 

Number  of  runs  were  tried  with  different  Initial  starting 
points.   The  approach  for  selecting  the  starting  points  was 
same  as  illustrated  In  the  previous  section.   The  details  of 
these  runs  are  presented  in  Table  26.   In  two  of  these  runs,  the 
value  of  the  projected  gradient  is  not  within  E^  but  the  value 


"  -■  9 
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of  the  objective  function  is  very  close  to  optimal.   A  typical 
convergence  rate  is  Illustrated  in  Sable  2?. 

The  optimal  solution  is  presented  in  Table  28.   The  solution 
to  the  same  problem  as  obtained,  in  [20]  is  shown  in  Table  29. 
The  value  of  the  objective  function  can  not  be  compared,  because 
different  expressions  were  used.   But  one  can  observe  the  closeness 
of  the  values  of  various  decision  variables  in  two  solutions. 

As  car,  be  seen,  the  gradient  projection  approach  proved  quite 
suitabl-  for  solution  of  this  problem. 

4.5   DISCUSSIOH 

Two  problems  solved  here  relate  to  the  situations  far  off 
from  those  encountered  in  practice.   But  the  approach  provides 
an  insight  to  the  problem  in  initial  exploratory  stages.   It  can 
be  seen  that  even  with  simplification,  a  mathematical  programming 
description  of  a  problem  can  retain  the  crucial  characteristics 
of  a  fairly  complicated  system. 

The  first  problem  introduced  an  approach  that  can  be  used 
when  one  or  two  tine  periods  can  be  isolated  from  the  rest  of  a 
project's  life.  This  method  was  then  extended  to  a  problem  with 
more  periods  and  involving  an  overyear  storage.  Both  these 
problems  involved  only  deterministic  aspects.  A  more  realistic 
representation  of  the  system  would  include  the  unpredictibily  of 
water  flows  and  other  random  factors. 

With  judicious  selection  of  various  limits,  tolerances  and 
initial  starting  points,  the  gradient  projection  method  proved 
quite  efficient.   The  convergence  was  quite  fast  and  the  solutions 


were  quite  accurate.   The  overall  approach  can  be 

superior  to  one  adopted  in  [20],  beoause  no  approximations  i 


CONCLUSIOH 

The  different  tost  systems  presented  in  this  work  suggest 
the-  usefulness  of  mathematical  programing  approach  in  the 
iv-  ■■■>'■  as;   and  management  of  water  resource  systems,  Ihese  exper- 
iments demonstrate  the  fact  that  the  crucial  characteristics  of  a 
fairly  complicated  system  can  be  retained  in  a  mathematleal- 
programming  description  without  redering  the  model  unduly  difficult. 

The  dissolved  oxygen  (DO)  model  presented  for  water  quality 
mo.na;;-3Si;rit  can  be  used  to  determine  the  minimum  total  cost 
associated  with  any  particular  set  of  minimum  allowable  DO 
concentrations  in  a  river  basin.   It  can  also  provide  the  useful 
information  regarding  sensitivity  of  both  the  cost  and  actual 
minimum  DO  concentrations  in  the  reaches  to  changes  in  minimum 
allowable  DO  concentrations  in  any  particular  reach. 

The  example  discussed  in  Chapter  4  were  connected  with 
water  quantity.   Though  the  problems  are  simple,  the  approach 
provides  an  insight  to  the  problem  of  water  resources  planning  in 
initial  exploratory  stages.   This  is  the  beginning  only.   If  more 
complexity  is  desired  the  benefits  from  the  flood  control, 
recreation,  urban  water  supply  etc.  can  be  incorporated. 

All  models  considered  were  deterministic  in  nature.   A  more 
realistic  approach  would  be  to  consider  the  stochastic  nature 
inhc-rant  in  inflows  and  water-demands. 

All  the  problems  fell  into  the  nonlinear  programming  class 
Characterized  by  linear  constraints.   The  gradient  projection 


method  developed  by  Rosen  [>]  proved  quite  efficient  in  solutio 
of  these  problems.   A  rapid  convereer.ee  rate  was  observed  in 
solutions  of  all  the  problems.  This  implies  computational 
efficiency  in  terms  of  computer  time  for  a  prescribed  accuracy. 

No  doubt,  the  success  of  the  method  is  largely  contingent 
with  a  judicious  selection  of  various  limits  and  tolerences. 
Also  the  results  obtained  are  very  dependent  on  the  design  of 
the  problem  to  be  solved.  Therefore,  it  is  very  essential  to 
have  basic  physical  knowledge  of  the  system.  The  program  can 
be  used  to  solve  problems  with  a  large  number  of  variables 
(about  60)  and  constraints  (about  150)  but  it  is  felt,  that  it 
will  be  more  efficient  for  fewer  variables  and  constraints. 
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APPENDIX  A 


GP  NOMENCLATURE  AND  SUBPROGRAMS 


NOMENCI, 

a-j  coefficient  of  input  constraint 

b1  right  hand  sta'e  of  a  constraint 

\*  •  t   of  the  inverse  matrix 

e  number  of  equalities 

F  value  of  objective  function 

S  gradient  in  the  direction  of  increasing  F 

H..  a  constraint  (hyperplane) 

3<  total  number  of  constraints 

m  nunber  of  variables 

HXNU  maximum  number  of  steps 

8XBN  maximum  number  of  re-inversions 

ni  constraint  vector 

Nk  constraint  matrix 

\  . 

(NTM  )~X  inverse   matrix 

q  q 

pS  projected,  gradient 

PD  projected  constraint  vector 

1  constraints  in  the  basis 

1  constraints  added  to  the  initial  basis 

u  linearly  dependent  constraints 

v  constraints  not  in  the  basis  with  X   =  o 

*  constraints  not  in  the  basis  with  X.  >    0 

x  variable  vector 

Z  unit  vector  in  the  direction  of  step 

P  gradient  interpolations  for  zTg  =  0 


gradient  interpolations  to  inert 

gradient  tolerance 

constraints  tolerance 

linear  dependence  tolerance 

interior  steps 

normal  distance  to  a  constraint 

step  counter 

temporary  flag 

step  length 


NAME  FUNCTION 

AKDA  Calculates  lambdas  {X  [::)). 

CLASS  Classifies  constraints  not  in  the  bas' 

and  W. 

COMMAT 
PUUCT 

MATCCX 


■i  :-:i,o.sf  mAnnrrs 


i 


Initializei      p,    y,   KOUiJT 
q,   q* 

Resdi     tJ,    t2,    6    ,    MXRN,    HXMJ 

na:c'    Ymax'7nax'  rma;c 

(7) 

Head i      Matrix  m,   k,    s,    e 

Constraints     + j ,    &1 . 

\U 

r?\ 

Headi      3     Vector 

-§■ 

X     Vector 

! 

I         Normalize   constrain1 


0 


INPUT  SECTION 


©- 

0- 


Add  H±   to 
all  i  in  e 

the  initial 

basis, 

i 

JT, 

mpute  all 
(x)  and  che 

U 

-© 


-© 


iHe-oheo^U)! ,@ 


<3>- 


.© 


-^"'ni^'"i 


© — B-^ 


3-Point  Extrapolatioi 


I    RSTUHF    |  I  List   all  1^    in  q*j 


O-S- 


.  11.   Flow  Chart  for  Subroutine  HEII 


To  test  HjL 

ENTSH.  1 »  adding  to 

basis, c'=l 


v  =  x,  k:v=x. 


KD 


MATRIX  COMPUTATIONS 
Fig.  12.   Flow  Chart  for  Subroutine  KATCOK 
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— 
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SsJ          "/ 
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•H     added 
NEXIT  =    3 


MATRIX  COMPUTATIONS 

Fie.    13-      Plow  Chart   for  Subroutine   CO:;:,AT 


APPENDIX   C 
COMPUTER  PROGRAM 


PROGRAM  FOR  THE  GRADIENI  PRGJECTII 

UN  X2UO),D<13,101,DN(10,10),A<25,1  >),GU  I,) 
P«1Uf.),n(25l,SD(2r.),-'P,(lJ),7  11O|,XC(10),Ga(lf 
XF(201,IH1 (10),R<10),AN8#*>25),Y( L0),PG( 1Q),1 

Y, PN, AfBOA , M,G ,X, F , P , EPS  1 1 , NEXlT.KQ, D, IH.NB.l 

r;PSI2,Ln,IU,R,KC;Q,MexiT,IV,KV,f,KMMi;,Il»,KW,K. 
,  I 'IT,  LIT. 

WATER/NFUNCKOUNI  ,K1 
!  ICI   (Y,PG),  (PN.PGNCRM) 

1=GRA0IE  a     tCLERENCE 

2  =  Cr.f:STKMNT  1CLERENCF 
3M.INEAR  DEPENDENCE  TOLERENCE 


05   FORI 

C'j2  FORI', 
CC3  FORI' 


f    OF  VARIAULES'6XI2) 
tf  OF  CONSTRAINTS  •  3X12) 
>l    PF  BOUNDS  '  8X12)  . 
«  OF  EQUALITIES'  5XI2I 
BOUNDS'  1 

COEFFICIENTS  rF  CONSTRAINTS' I 
NORMALISED  COEFFICIENTS  OF  CONSTRAINT! 
NORMALISED  LIMITING  VALUES  OF  CONSTRA! 
LIMITING  VALUES  OF  CONSTRAINTS') 
INITIAL  X-VECTOR' I 


I  FOR  NONLINEAf 


AND  PRINT  LIMITS  AND  TOLEREf 

12,EPSII,FPSI2.EPSI3,TMAX 

AX6,MXRNtMXf 

,cl>SI2,i 


KOUNT-C 


K1  =  0 


HUN-.C  .  !1*THAX 

•  NEK  MATRIX  '    READ  C, 

RE  ".i:  1    ,M,K,NB,NI 
PR  If 


N8P1=I 
KM  *B=I 

ifinb: 


,24,22 

AND  PRINT  SUBSCRIPTS  FOR  BOUNBS 

%  IJXB(  I  ),I=l,Ne I 


REAC  AND  \C^KALIJE  CGMSTRA  I '-IT  S 

,D  K'iSIHI  l,J),  J=1,M),I«1,KHNB) 

NT  10  ".7 

^T  1040,  <(A(I,J),J=1,M),I=1,KPNB) 


SIM  I  )=;(.:T(;,I'U  I  ) 
DO  26  J=  I  ,'■: 
>  Al  I,  J!--A(  I  ,J)/SOI  : 


ACD  II  TO  INITIAL  RASIS  FOR  ALL 


2  7  If  CIE)  30, 30 1 28 
28  K«=0 

IH*NB*1 

DC  61C  l«l,NE 


640  PN=PN4Y(  I)  :■  ,' 

PN=SC«T(PN) 

IF(PN-CELTA)  650,650»-64! 
645  DELTA=PN 

INDEX-N 
650  IFIN-NEI  625,655,655 
655  IFUNCEX1  675,675,660 
66   IH=INCEX+NI 

CALL  CATCOI  !1  I 


DO  6t'C  1  =  1, KC 

DO  680  J=l,KC 

680  DM  « , J ) * D ( I , J I 

RFAC  P-VECTO* 

1   READ  1040,  fB(I),I=l,Kl 

PRINT  1040, (BUI,  1  =  1, K) 

IF(KfNE)  4:<,^:,31 
31  00  32  I*1»KMNB 

3c1  3(J)=0(J)/SD(I) 

PRINT  1009 

PiUMT  K.40,(P.(I  1,1  =  1, K) 
4i  CONTINUE 

41  NFUNC=: 


COf'i 


A(X 


52  CALL  AfOA 

»F(NB]  54,54,64 

64  DO  5  3  I  =  1,W 
IF<AMBCA(I»*E>S12)  59,53,53 

5J  CONTINUE 

54  IF(NBP1-K>65,65,130 

65  DO  57  J*NBP1,K 
IFU-NBP1-NE)  55,56,56 

55  IF(ABS(AM8DAU)  I-EPSI21  57,57,59 

56  IF(Af6CA(J)+EPSI21  59,57,57 

57  CONTINUE 
GO  TC  13-' 

59  PRINT  1040,  IXU),  I»l,M1 

Ml^f  1040, (AfBDAIJ), J=1,K) 
IF<KC)  80,80,60 

6(  DO  6  1  I=1,KC 


CCS       CTIGN         COMPUTT 


7.    DO    71     J=1,M 
71    VU)  =  XU) 

CALL     MiCO!    (  i 


-El'SI2)     62,62,110 


1HK-2)    89,8, 
86    DO    88    1*2, K 
IFIAMBCAI  I  I- 


90    CALL    PATC0M1) 

GO    TC     HOO.l-CJSJiNEXIT 
95    KQMl=KC-l 

IF    (KCM)    97,97,94 
9'.    CO    96    I»1,KQM1 


:    f F ( S I G i» A 1    !!',r.i,n: 

L  C\LL  CCMMA1  IEPSI3I 
GO  TC  1105,1031 .MEXIT 


NC    FEASIBLE  ) 
MSI  5001 


NOT  FEASIBLE  ) 


RE-CHECK    LAMBD) 


•RINT    : 


,(X(J 


NT    5.V 


5004    FORMAT     (•       LAMMS'  ) 

P.UNT     1040,  (Af  TDA(  J)  ,J=1, 
CALL    FUNCT<X,F,G,KQ) 
IF(KGLM)     133,133,140 

133  IFIKC-KFQ)     9999,138,134 

134  CONTINUE 
130   CONTINUE 


ETA 


.    CLASS 


CLASSIFY 


FOR 


G»#ORW»C 
DO    14  1    J*1,M 
141    GNQRM»GN0RM+GIJ)**2 
GNORM*SORTlGNORMJ 
PRINT    5006tRNORM, IG(J1,J»1,M 
5CC6    FORMAT!'  GRAO  I  EM  ■  ,  F  12  .  6/ 

IFIGNCRM-EPSU1     142,142,150 


C  PRCJOCT    GRADIENT 

C 

17      DO    169    J*1,M 

164    V(J)=G(J) 

call   catco?  (2) 

PGMORC*a.C 

DO    171    J=1,M 
171    P  GNftR ^ ■ PGNO RM+PG ( J ) *  *2 
PSNORF=SORT(PGNORMI 
PRIN1    5007.PGNORM, IPGU),  J=! 
5C07   FORMATS      PP.rj   GRAO'  .F12.6/I 
IFIPGNCRW-FPSU )     175,175,17; 


1,173,180 

:  HOT    ZERO  .\T 

I  NUT  ZERO  AT  \ 


IF(i 

EXI 

' 

745, 

745, 

4  60 

ECL 

ss 

r-VH 

FOR 

t  MOT 

IN  Q 

GO  T 
75  IF(S 

C  1 

? 

76,177,17 

6 

P 

KCJ 

CT 

ON  It 

RO  A 

FTER 

ORCP 

09  FORf 
CO  1 

77  IF(K 

78  CALL 

GO  T 

T  5 

AT 
C  1 

CO 

0  ( 

Q) 

20 

PRCJ 

9999 
(EPS 
179) 

ZERO 
420, 

AFTE 
178 
IT 

R  DRO 

IFIKECP-KC)  1903,1903,1902 
1903  CO  19C1  J-KEOP.KQ 
1901  RRU)«RU) 


914 


GO  TC  1925 
l»ZU 


GO  TO  19?5 
1915  KK=IVU)-NB 
ZI»~Q.C 
CO  192  J  =  1,M 

192  ZN=ZN+Z{J)*A(KK,J1 

19,?-.  IF(ZN-CELTA)  19  3,1935,1935 

193  INDEX* IV (I > 


97  IF(KECF-KQ1  19" 


946 

LLL=  : 

GO  TC 

23T 

200 

[F(N: 

TA)  2C1, 

20  A 

201 

201 

[FlKf, 

-KE0I20; 

,?v 

t2C 

202 

CALL 
GO  H 

CCVMATIS 
l2r>4,2C 

3)f 

IEXI 

203 

NETA  = 

C 

FORMAT KX2H    Q,I2,3H    F 
IFIKC>212,21?,211 
PRINT    9211,111-1  U),J' 
1F(LUC)214, 21^,213 
PRINT    9213, (IU( J),J=1 
IF(KV)216, 216, 21=5 


[FIJ-M>)232,232,235 


!F(  n-[)238,23' 


PRIN1  5G11.K0UN1 
501 1  F03MA1  (•  STEP"  1 

140  on  ?'.  1  ,i  =  i,;' 

241  X(J)=Y(J)+l*ZtJ] 

KMU0=PGNORN 
2-..  C.Ll  FL'lCTIX.F.C 


,KU 


00  251  J=1,M 

251  ZG=ZG  +  ZU)*G(J) 

P;U-;T  5;i2,MU,IMT,NETA,L,T,F,ZG 
5012  F0RMA1  <4X5l'"cI.".  =  ,3XI3,4X6HGirNA=,3Xl3,4X5H''iFTA=,3XI3,3x;: 
lX2HT=,F12.6,4XtftF  =  ,F12.6,4X3HZG*,F1.2.61 
IF(T-Tf-IM)270,252,252 

252  fF(MU)S999,253,257 

2  5  3  [F (MUPAX) 9999, 27C, 254 
25',  IF(  JG+EPSI]  )-?55,255,270 

255  DO  2f,6  J  =  l," 

256  X2(  J)=Y(.J) 
GO  TO  B255 

2  57  I F ( HU-f UH*X 1258 ,270, 270 

258  tF(ABS(ZG)-EPS! 1)8255,8255,8250 


IF(ABS(AM8DA(J) I -6P5I 2)281 ,281,300 

L  CONTINUE 

,  DO  2  6  2  1  =  1,1- 

T  F  (  A,-  BCA(  I )  +EPSI2  1292,282  ,282 


Riio-Apatmi) 

IFIK-21298, 2^9,299 
)  DO  295  T=2,K 

[ F I  AMBCAd  )-RHO  1 294, 295 , 29 


CALL    KATC0M1I 

CU    10(309,300, 296), NEXIT 

296 

NETA*C 
INV*1 

30  C 

60   TO    310 
IFCSIGKA)301,31  ),301 

5021 

FORMAT!'         CONSTRAINT    VII 

31  : 

siV.ma  =  siv-:a  +  i.O 

5C32 

FORMAT!'          X-CORRECTION' 

CO    311    J«1,M 

311 

CALL    CATCQHt3) 

312 

DO    312    J=1,H 

X  !  J  )  =  V  (  J  ) 

INV-1 

r.O    TC    290 

320 

IF(SIGPA)321,325,321 

321 

CALL    FL'ICT(X,F,G,KI,') 

IF(F-FY)322, 325,325 

5013 

FORMAT! ■            F      DECREASE* 

lF!NEXIT)325,325,46Ci 

3  2(: 

PRINT    5014 

51  ]'<    i  0  IMATI'         / ^ A X I 

go  ro  /if 

327    PRIN1     5    15,<X<J 
■     15    FORMAT I6F12. 6) 

GU    TU    IV 


riON    SCHEME 


(J)=LCI.X 
G=7G4Z( J 

3n=zcc+z 

t  (ZGCI3 
F<ZGC-Zi 


\)332,332,33- 


J)*Cri(J) 
,33b, 338 
33"5, 335,336 


341     K2<J)*X(J) 

PRINT    'Cl^.r.ZG, !Z(J),„=1,") 
5C16    FU''MMC  3    PT.     INTrRPCLAT ICV 5X2hT=F 

1    5X3HZG»F12.3/'  Z-VECTOR' / ( 6F  12.6 ) 

GU    TC     351 
35  ■    r»TPRIME 

PRINT    5.  17,  r,7G,(Z{j)  tJ'liM.) 
5017    FORMATC       3    PT.     EXTRAPOLATION' t5X2HT=F 
1    5X3HZG«F12.3/<  Z-VECirv  /  (6F12.6  I 

351  IF{T-Tf'AXI353,353,352 

352  f=TMAX 


'ECTOR    GRADIEt 


5020  FORMAT (  13 


' 


DIMENSION    X  ?  t  1  ?  )  ,  D  (  1  1 , 1  ■'  ■ )  ,  r ,  J  (  1j  ,  1 0  )  ,  A  (  2  5  ,  I  o  ) ,  G  1  : 

I  IU(1C),1W(?5)  ,l'<25)  ,SD(2r.)  ,RR(1  J  I  ,  ,'(10),XO(10),( 

I I  (  1 "  )  ,  .•  X."  (  2  "  )  ,  I  H  I  (  10  )  ,  P.  t  10  )  .  AMtSDA  ( 25 ) , Y ( 10  > , PG ( 
IBEtXU    ) 

CQWHCIS    VfPN.APBOA.M.G.X.FtP.EPSIUNEXIT.KQ.CIH, 
"  2tLD,IUiR,KE0tMEXITiIV,KV,B.KMNB,IW|l 


1  3,23 


3': 

1F(KC(. 

>  36,  36,31 

DO    3  5 

1=1, KEC 

DU    3  5 

J=1,KI=Q 

36 

n)l'cJ) 

KO«KEC 

PRINT    935 
935    FORMAT (•  NEW    BASH 


IF     (KbC)     54,54, 52 
:    DO    53    1*1, KEG 

DO    53    J=l,KFC 


INT  =  0 

LDC  =  C 

Ka=KfcC 

PRINT    935 

INV  =  C 
55    CONTIMiE 
60    IH=IIJI(K"*1) 

CALL    MTCQI  i  1 


DIMENSION  X2(1J),D(1C,10) iDNI 
V<  LOI  ,JXB{2r  IrlHI  liC  )  ,R(IC  i  ,A 
IU(lC),IWt25J,e(25>,SD(25),RR 


>SI2,ld,iu,r,keq,> 

LMU,NETA,JNT,LDC 

FQlll  «.M-MC:t      (Y,I-'C.)  ,  (PN.PC 


;    CCMP'JTAf IONS    I 


NEXI1 .KQtDi IHiNBtV,A,JX6i [HI, 


PN=l.C 
GO    TO    261 
3^     1HIH-N6) 


i  «0  50   I  •    lil 


GO    TO 

TC   vun 

GO  TO 
8C  I  FLAG 
90    IF     IK' 


DO  150  I  =  l.KQ 
KK  -  II  It  I 
IF  (KK)  15 


co  rc  150 

1  IF  IJ-JK) 

I  Y(J)  -  Y(. 

GO  TC,  150 


li(L)  LV 
PRINT  9200, 
1  FORMAT!'     H  ',\2,'    LINEARLY  DEP1 
RETURN 

:  j  =   kc  +  l 

DC1.J)  =  1.0/YB2 

IF  (I-JI  24C.240.230 

1  =  1-1 

D(I,J)  •  D(  I.J)+RU)*R(  I  1/YE2 

D  (  J  ,  I  )  =  D  (  I  ,  J  ) 

GO  TC  22' 


IHIIKCI  =  IH 
H(L)  ADDEl 
"  9260.  1 


926  -  FORMAT ( 


ADDED    PN«,F12.6I 


SUBROU  IME  CCI 
THIS  SUPROI 


136 

RY  MATRIX  CCHPUTAT1 

S 

25, 1C >,GllO>,XtlO)i 

PUO), 

),XOUC),GOI10>,Xl( 

VI  10),PGI  LOW  IV  (101 

r,KQ,D,IH,NB,V,A,JX 

B.IHl, 
INV,ON, 

42  DELTA 
«  IF  (L- 

KOI  « 

GO  TC 
45  IF  INZ 
«6  MEXIT 

RETURfv 

)  46 

50  CO  60 

I  = 

IH  =  I 

Mlli 

NiZ-KO)  63,63.1 


(I)  =  IHI ( If  11 
(n;c)  75,75,65 
7o  J  =  l.NZM 

: ,  j  i  =  D  ( I + 1 ,  J ) 


105  DO  13C  I  ■■ 
DO  13C  J  : 
IF  (I-J) 


12.  0 ( I , J  )  - 

13.'  CONTINUE 

IF  (ICC) 


SUBROUTINE  AMOA  1 

THIS  SUbROUTINE  CARRIES  OUT  CALCULATION  OF  LAMBDAS 
DIMENSION  X2<10),D(10,lu),DNtlO,10),A(25,10),G(10),X(10),PUO), 
11 'J  I  1  ),H  (;■■"■)  .'  (?'.  l.SUU'*  )  ,  -I'M  I  -  !  ,  /  !  i.  i.  *Oi  1  ;)  ,GU(  1  •.■')  .XI  (1  J)  . 
IVI 10) , JX8I2C I ,IHI < 1C ) ,R( 10) , AMBDAC25) , Y 1 10 ) .PG11G ) , IV ( 10  >, 

"r   |KHCN    y.PN.AMBDA.M.G.X.F.P.EPSII.NEXU.KQ.D.IH.NB.V.A.JXB.IHI, 

I    EPSI3,LeSI2,LD,tU,R,KE(M-'EXlT,lV,KV,B,KMr.,U',KW,K,MXiW,l\V,r    . 
lMU,NETA,tNT,LCC 
[Fir  i     ,EQ.     r  )G0    TO    '.0 
10    DO    3C     1  =  1, NB 
J=JXBU) 

1FU   .ir.   osn  to  i'C 

AHBOAd  l«X(J)-B(  I) 

CO    TO    30 
20    AMBUAI [)=-X(-J)-B< I) 
30    C0M1  IMC 
4  I    IV  IKI  NB)8(  ,80,5  ' 
50    00    70    I=1,KMNB 

TOTA--C. 

KK*NB+I 

DO    6C    J  =  1  ,  f '  ' 
67    TCITA  =  TCTA  +  A(  I,JI*X( J) 
K     AMBDA(KK)*TOTA-BIKK) 
Bl   RETURN 


TIIS     SUK<Jl..TI\i;    CLASiiriES    THE    Ct\'SrKA 
DIMENSII    I    X2U0)  ,;;(!•  ,1  .)  ,DN4 1SV10 ) , A (  ■>:,, 

1 1 ui  i c, ) ,  I .:  ( z'-> )  ,t  ( .? ■> ) ,  sr.  ( ;■•-. ) , rr  i  i ) ) , i  ( i o  , ,  X; 

1VI10)  ,  IXB(2<  I.IHI  (10)  ,R  (  H    I, A  II   .'.(      >1     Ml 


XIT.KO.D, IH,N8,V 


DU    20    J  =  l ,KU 


AIER    CUAL 1 1 v     ' 

CULATES'TI  ;   VALl 


1 

D 1 1 

10       FORMAT ( 

IC1    1  '•'  t>  •'  1  ( 
!     2    FORI  AT 

N    X  ( 1 '  1  .  G  !  1 

[1         11 

1  1  2  .  '  /  ' 
11'!,'            CO 

□i    r,Ki 

1  ,  C  1  (  i  ■•■ )  ,  C  2  1  1 

I.  •■    •    *    '  1 1  i .  ' 
nstraim:    in 

-VALUE' ,F12.6 
5T    COEFFICIENT 

i  .  .       ■  ■ 

«    OF    FONC 
BASIS    =     •  ,  H./ 

S    Fn;<    PLANTS'  ) 

1,01(101,0 
IONAI    EVAt 

**CCST,i 

-  C    AND    EC     INFO 

IMATION**' ) 

LC«    FORM/   I 

•       REACH    S 

1    :  ;■  !  1  •' 1   M 

COST                 K 

HX.BOC         Ml 

.  .,    FORMAT  (7X1 2, 5F1 2. 6) 

FORMAT     (•  MIMtKUf    TCIAL    CCST'  ,F1J 

iCOC    FORMAT  (3F15. 71 

.'  A     II    ••■AI(I'i,3F15.7» 
0RMAT(7F8j3) 

IFlNFtNCJl ',!',..  ' 

REAC    AND    PRINT    CCST    COEFFICIENTS 

1       READ    1CI  ',  (CHI  t,C2tU,C3(l  ),  I«1,K1 

P  R [NT    1    2 
DC    11     1*1 tM 


CHI  I  >*Xtf»*2. 

LM.'.FL'X.KC 
(I».lxl,KI,F 


JRADIENT     VECT< 


rfac   r nr:  and   co   ini 


2A    B2UI 

CA 

t:  1  t  1  ) 

•  U.-X3I  I)  )    r  i  |  I  ) 
LCULATt    &C  rUAL    CO    11    P 

FACHC 

D 1  ( A  ) 
DK5I 
114) 

0  11.'! 

dii  n 

I  *  X 1 4  1 

.  I     j  151  "     1  ' 

•-1.9P8039  +  C.  154961  *X< 
i.324243*Xf51 

16.7] 
U4.5S  ?280*XJ5)-t-2.895C 

U  .  5  3 
1  +  1.2 

75564  U  2 1 
526  )23*X( 

t  crsrs 

)**2 

.. 163363* 
.155779* 

+  .2',  32  3; 
N. 41732 

X(i 

CO    27 

n  ciin 

)  =  CIU)*Q2(I  )*X<  I  HC3I 

I  )*X( 

CO    2tf 

28  cm  i 

<]  )  =  C1  (  I  )*C2(  1  K-x(  I  HC 

3(1)  * 

( I  )**2 

3)+4.35!  409 


CALCULATION    OF    TCTAl    CCS 


oltput  cnsT,(-rc,nr   information 


IC5,I,X3( I),C11(!),C2(1) 


SUBROUTINE  FUNC1  (X,F, 

TWC  PERIOD  PROBLEM  FOR  WATER  P.r 


THIS  SUBROUTINE  CALCULATES  THE  FUNCT 

ION  VALUE  Ar 

DIMENSION  X(2C).G(20) 

C0MMCN/WATER/NFUNC,K0UNT,K1 

100 

FORMAT  (•      ITERATION  //  =  ',13,'                # 
1=  ',I3t'      #  OF  CONSTRAINTS  IN  BASIS 

OF  FUNCT 10 

101 

FO  1ATI5F12.6/'     F-VALUF'  ,  F12.6  ) 

lCGl 

FORMATdH-,'     BENEFIT  FROM  ENERGY 

•,F12.6) 

1C02 

FORMATdH-,'     BENEFIT  FROM  IRRIGATION 

' ,F12.6) 

1003 

FORMATdH-,'     COST  OF  RESERVOIR  B 

•,F12.6) 

1004 

FORMATdH-,'     COST  OF  RESERVUI'.  C 

•  ,  F  1  2  .  6  ) 

10  1C 

FORMATdH-,'     MAXIMUM  NET  BENEFITS 

IFIK1)10,10,20 

NFUNC»NFUNC+1 

•,F12.6) 

CALCULATE  THE  FUNCTIG 


ILUE  AND  THE  GRA0IEN1 


453. 7*AL0Gd. +0.2*1 


/(l 


Gdl=2.*Xd>*229.4 

G(2)  =  l. 4+1453. 7*0. 2/d.+0.2*(X(2)  +  XC 

GI3)=-48.7+( 453. 7*0.2/ d.+O. 2*1 X(2)+) 

G(4)=-43.*(l./(l.+.2*XI4)))+43.*(.2*) 

G(5)=-47.*d./ll.  +  .3*X<5) ))+47.*I.3*> 

PRINT  100,KOUNT,NFUNC,KQ 

PRINT  101,  tXU)  ,I  =  1,M),F 

GO  TO  40 


,TE  BENEFIT 


AND  COST  FUNCTIONS 


X(8)=43 


20  X(6)=229.4*X(l)+Xtl)**2 

X(7)d.4*X(2)-48.7*X<3)+453.7*AL0Cd.+0.2*(X(2H 
X!4)/d.  +  .2*X(4> ) 
X(5>/< l.+.3*X(5) ) 

PRINT  1001, X(6) 

PRINT  1002, X(7) 

PRINT  1003, X(8) 

PRINT  1004, X(9» 

PRINT  1G1C.F 
40  RETURN 

END 


:  iji  Ci  IX,  F  ,G,KG!  1*3 

FOUR  PERIOD  PROBLEM  FOR  HAt6«  RESOURCES 

THIS  SUBROUTINE  CALCULATE?  Flit  FUNCTION  VALUE  ANO  THE  GRADIENT 

DIMENSION  X(20),G(20),XXI(20)tFH10»,F2(10l 

CilMr'CH/l.'Al  ER/NFUNC,KUUNT,K1 
.  FORMAT  I'      ITERATION  K    =■    ',13,'      »  OF  FUNCTIONAL  EVALUATIONS 


02  FORMAT! IH-,  •  PERIODS     1        2 

03  FORMATUII-,'  STORAGES  ',  4F8.  3  ) 

04  fli :'.:  :.\T  (11!-,  ■  OUTFLOWS', ',F8. 3) 

05  FORMAT! 1H-,'  ENERGY • , 4F 8. 3) 

06  FORMATUII--,'  TOTAL  FNERGY  •  ,  F8  .  3  ) 
!C7  FURMATUH-,'  CAPACITY  OF  RESERVOIR' 
108  FORMAT  !  It!--,  '  BENEFIT  FROM  ENERGY', f 
:09  FORMATUII-, •  COST  OF  RESERVO  IR  '  ,  F8. 
!10  FORMATUH-,'  MAXIMUM  NET  BENEFITS', 


CALCULATE  THE  FUNCTION  VALUE  AND  THE  GRAOIENT  VECTOR 
F--2.2SA*X(6  )  +  .01*<X<6  )  **2  )  -.43*X  I  5  )  /  1 1  .  +  0  .2*X  (  5  )  ) 


CC.  1=0.0 

G(5)=-.A3*(1./(1.+.2*XI5) ) 1+.43* 
G(6)  =  2.29<t+.02*X(6> 
PRINT  100,K0UNT,NFUMC,KQ 
PRINT  101, (X(I) ,1=1, M) ,F 
GO  TO  40 

READ  INFLOWS  TO  THE  RESERVOIR 


,TE  OUTFLOWS  FROM  THE  RESERVOIR 


DO  22  1=1,3 
!  XIII  I*A»  =  X1U  ■ 

X11(8)=X11(4) 


REAC  ENERGY  I 
READ  1001, IF2U 


CALULATE  BE 

.  l-  i  ;  r 

H16|-.':3*X1 
INT  1002 

LI 9 1  I 

INT  1003, (XI 

till 

INT  1004t(Xl 

'('■). 

INT  1005,1X1 

INT  10Q6,Xll 

(10) 

INT  1007, Xll 

UNT  1008, Xll 

(  15) 

INT  1009, Xll 

(  16) 

INT  1010, F 

iMD  cusr   I 


l.+0.2*Xll (91 1 


OPTTMI'.'.ATIO" 

OF' 

\:\:vi-.  ].■'  '.en, \  \ 


PRATAP  BHANJI  DAUH 


AN  ABSTRACT  OF  A  KASri^.'S  THFSIS 


submitted  in  partial  fulfillment  of  the 


requirements  for  the  -eartf 


MASTER  OF  SCI?  :CK 


Deaprtment  of  Industrial  Engineer!  rs 


In  recent  years  at  tempts  have  been  made  to  apply  the 
various  mathematical  programing  techniques  to  the  planning, 
design  and  operation  of  water  resource  systems.   Two  general 
types  of  models  have  been  fruitful  in  the  field  of  water  resource 
development!  the  simulation  model  and  analytic  model.   In  this 
thesis,  two  deterministic  models-one  connected  with  water 
quality  and  other  connected  with  water  quantity-are  proposed 
and  solved  within  the  framework  of  an  analytical  approach. 

For  each  system  that  is  considered,  the  procedure  involves 
the  understanding  of  the  basic  physical  systems;  the  development 
of  systems  equations  or  mathematical  models  and  the  solution  of 
the  i'V  posed  models.  Data  used  in  various  models  was  drawn  from 
the  literature  wherever  possible.  This  practice  was  adopted  to 
insure  the  realistic  response  behavior. 

Like  most  of  the  models  used  to  describe  the  water 
resource  systems,  the  models  presented  here  are  characterized 
by  a  nonlinear  objective  function  and  linear  constraints. 

Rosen's  gradient  projection  method  appears  to  be  a  powerful 
tool  for  this  class  of  nonlinear  programming  problems  characterised 
by  linear  constraints.   The  specific  purpose  of  this  thesis  is 
to  apply  this  technique  to  the  various  water  resource  models 
and  analyze  the  results  to  derive  optimal  design  and  operation 
policies. 

The  method  is  initially  described  in  some  detail  with  an 
eraphasis  on  conpv.fcaticnal  aspects.   Then  the  method  is  success- 
fully applied  to  solve  the  various  models  describing  water 


The  construction  ,'  id  sclu!  i  >n  ol    r  J   1   ati  caJ  tnc  r]  oj" 
dissolved  oxygen  for  a  8imp].ifled  rive  5  basin  indicates  how 
a  mathematical  model  can  generate  useful  water  quality  control 
information. 

Next  a  simplified  river  basin  configuration  is  selected  to 
illustrate  the  mathematical  programming  approach  to  water  quantity 
aspects.   Such  an  approach  is  useful  in  initial  exploratory 
stages  of  water  resources  planning. 


